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EXECUTIVE SUMMARY 


The Lower West Coast (LWC) Planning Area of the South Florida Water Management District (SFWMD) 
faces numerous water management challenges. Growing freshwater demands, dwindling traditional water 
sources, increased levels of environmental protection, and sea level rise need to be addressed and managed 
to protect the area's water resources and provide adequate water supply. Regional water supply plans are 
the SFWMD’s primary tools to address these issues. 


To evaluate potential impacts to water supply and provide a tool to support water supply planning, a three- 
dimensional groundwater flow model of the surficial and intermediate aquifer systems underlying the LWC 
Planning Area was developed. The Lower West Coast Surficial and Intermediate Aquifer Systems Model 
(LWCSIM) was designed to simulate the response of the aquifers and natural systems (e.g., wetlands) to 
stresses such as proposed wellfield pumping. Results of the model applications can provide guidance in 
developing water management strategies, support periodic updates to the regional water supply plans, and 
be used in regulatory applications. 


The LWCSIM was developed using the United States Geological Survey's MODFLOW model code. The 
geographic extent of the model covers extends from Charlotte County down the Gulf Coast to Everglades 
National Park in Monroe County, and from Lake Okeechobee in Glades County to the northwestern corner 
of Miami-Dade County. This area was divided into a uniform rectangular grid, oriented north to south, with 
grid cell spacing of 1,000 feet. The model comprises nine primary layers, including five productive aquifers: 
Water Table aquifer, Lower Tamiami aquifer, Sandstone aquifer-carbonate zone, Sandstone aquifer-clastic 
zone, and Mid-Hawthorn aquifer, and four semi-confining zones between the aquifers. 


The model was calibrated to steady-state and transient conditions. The steady-state model was calibrated to 
average 1999 water levels, which ensured the water budget was reasonable. The steady-state model 
represents the approximate conditions of 1999, which was a period of minimal change in storage in the 
study area, and a period of increased monitoring that provided the initial conditions for the transient model. 
The transient model was calibrated to the period from January 1999 through December 2012, a total of 192 
monthly stress periods, with a validation period from January 2013 through December 2014. 
Approximately 500 measured groundwater levels were used as head calibration targets. Four measured 
structure flows and estimated baseflows were used as flux calibration targets. Manual and automatic 
calibration (Parameter Estimation software [PEST]) methods were used iteratively during the model 
calibration process. The primary purpose of the transient model is to provide a tool to conduct long-term 
(20 to 50 years) planning and evaluate water supply issues. The LWCSIM is designed to evaluate regional 
conditions from stresses on the surficial and intermediate aquifer systems in the LWC Planning Area. 
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Next Generation Weather Radar 
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Natural Resources Conservation Services 

per capita use rate 
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Pearson coefficient of correlation 

surficial aquifer system 

South Florida Water Management District 
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Tamiami confining unit 
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Upper Floridan aquifer 

United States Army Corps of Engineers 
United States Geological Survey 

Web-based Hydrologic Analytical Tool 
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1 INTRODUCTION 


This document presents the Lower West Coast Surficial and Intermediate Aquifer Systems Model 
(LWCSIM), including its purpose and development. The LWCSIM boundary (Figure 1-1) encompasses 
the Lower West Coast (LWC) Planning Area, one of five water supply planning regions in the South Florida 
Water Management District (SFWMD or District). Тһе LWCSIM was created to assist with long-term 
(20 to 50 years) planning and evaluation of water supply issues in the LWC Planning Area. The model may 
also be used for regulatory applications, which could require modifications to the model to meet regulatory 
criteria. The LWCSIM was developed following standard industry protocols (Franke et al. 1987, Anderson 
and Woessner 1992, Reilly 2001, Reilly and Harbaugh 2004). 
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Figure 1-1. | Model location and boundary. 


1.1 Water Supply Geography and Issues 


The LWC Planning Area includes all of Lee County, most of Collier County, and portions of Hendry, 
Glades, Charlotte, and mainland Monroe counties (Figure 1-2). The region generally reflects the drainage 
patterns of the Caloosahatchee, Imperial, Estero, and Cocohatchee river basins. The LWC Planning Area 
includes many coastal and inland natural systems, including Big Cypress Swamp and a large part of 
Everglades National Park. 


The LWC Planning Area faces numerous challenges to maintaining adequate water supply for growing 
urban and agricultural demands while simultaneously meeting the needs of the environment. Detailed 
information on estimated and projected demands for the LWC Planning Area can be found in the most 
recent update to the water supply plan (SFWMD 2017). 


Historically, water demands in the LWC Planning Area were met using surface water and groundwater. 
Surface water has been the primary source of water supply for the agricultural industry in the 
Caloosahatchee River (C-43 Canal) watershed. However, surface water from the existing canal and storage 
networks alone is insufficient to meet agricultural water use demands and environmental needs during 
1-in-10-year drought conditions. Past analyses concluded that additional storage was necessary to provide 
adequate resources to meet existing legal user and natural system needs in the LWC Planning Area. 


Throughout the LWC Planning Area, the surficial aquifer system (SAS) historically served as the major 
source of fresh groundwater for Public Supply (PS), Landscape/Recreational (L/R), and Agriculture (AG). 
However, analyses of the SAS indicate it is a limited water source in many areas. 


Saltwater intrusion is an ongoing concern, resulting from continued use of shallow groundwater sources 
near the coast, potential sea level rise, movement of connate water in the SAS, and upconing of poor-quality 
water from deeper aquifers. 


Historically, the Sandstone aquifer (SSA) and Mid-Hawthorn aquifer (MHA) within the intermediate 
aquifer system (IAS) have been important sources of fresh water for portions of Lee and Hendry counties. 
However, the freshness and productivity of these local aquifers are not consistent throughout the LWC 
Planning Area. Analyses indicate these localized aquifers are limited water sources in portions of the 
planning area due to the cumulative effects of withdrawals by all water users, which decrease water levels 
in the IAS and could harm the resource (e.g., through saltwater intrusion). 


To reduce the potential impacts to traditional freshwater sources, future increased demands are projected to 
gradually shift to alternative water sources, including increased use of the underlying Floridan aquifer 
system and reclaimed water. The Floridan aquifer system is a regional resource, underlying portions of 
Georgia, South Carolina, and Alabama as well as the entire Florida peninsula. 


1.2 Objectives 


The objective of this work was to develop a calibrated groundwater flow model that incorporates the SAS 
and IAS in the LWC Planning Area to support the regional water supply plan updates. The intended 
purposes of the model include the following: 


1) Facilitate the SFWMD’s evaluation of the current and future availability of groundwater in the SAS 
and IAS; 

2) Assess future water supply and management strategies; and 

3) Evaluate potential impacts of groundwater withdrawals on wetlands and natural areas in the LWC 
Planning Area. 
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Figure 1-2. South Florida Water Management District planning areas. 


1.3 History of Modeling in the LWC Planning Area 


Тһе SFWMD retained Marco Water Engineering, Inc. (2006) in April 2005 to develop a MODFLOW-based 
groundwater flow model of the SAS in the LWC Planning Area. The Marco Water Engineering, Inc. (2006) 
model, covering 5,129 square miles, simulated the Water Table aquifer (WTA), Lower Tamiami aquifer 
(LTA), and SSA. The model grid consisted of 765 rows and 622 columns, and the cells were 704 feet (ft) 
by 704 ft. A transient model calibration was conducted for a 15-year period (January 1986 through 
December 2000). Table 1-1 illustrates how the LWCSIM layers relate to the Marco Water Engineering, 
Inc. (2006) model as well as nearby models (Lower East Coast Subregional Model [LECsR; Giddings et al. 
2006] and Districtwide Regulation Model [DWRM; Environmental Simulations, Inc. 20141). 


Table 1-1. Model layering nomenclature across multiple models. 


Marco Water 
Aquifer ee LWCSIM Model рше ү. Lower East Coast Subregional 
System Mode Layer Name в шм Model Layering 
Layer # Model Layer 
Number 
Surficial ; wae able Тілі | Layers 1 and 2: Holocene + Q5 
amiami confining 
Aquifer 2 unit Water Table EOE OST Oe OT 
Syst Tc 
DE 3 Pov ш 2 Layer 3: T2 + TI 
aquifer 
4 Upper Hawthorn 
confining unit 
5 SSA clastic zone 
| 6 SSA confining Peace River 
Intermediate layer Aani 3 
: quifer (PZ1) 
Aquifer 7 SSA carbonate 
System zone 
8 Mid-Hawthorn 
confining unit 
9 Mid-Hawthorn | Upper Arcadia 
aquifer Aquifer (PZ2) 


LWCSIM = Lower West Coast Surficial and Intermediate Aquifer Systems Model; SSA = Sandstone aquifer. 


The LWCSIM was developed to incorporate a) the IAS in addition to the SAS, b) 15 years of additional 
water use and water level data, and c) updated hydrogeologic data and hydrostratigraphic interpretation of 
the data. It is important to periodically update models with current data. The improved model will facilitate 
evaluation of the current and future availability of groundwater in the SAS and IAS. 


1.4 Peer Review 


An independent, scientific peer-review panel evaluated each major phase of LWCSIM development and 
calibration. The following groundwater modeling experts with experience in Florida composed the 
peer-review panel: 


е David Е. MacIntyre, Р.Е. (Chair), AquaSciTech Consulting, PLLC 
e James O. Rumbaugh, P.G., President, Environmental Simulations, Inc. 
e Mark A. Ross, Ph.D., P.E., President, Hydrosystems, Inc. 


The peer-review panel was engaged from the conceptual model development phase through calibration and 
model documentation. The first panel meeting occurred in October 2016. Throughout the peer-review 
process, publicly noticed teleconferences were conducted periodically to update the panel on the LWCSIM 
progress and solicit input via a series of tasks outlined in the Scope of Work. All communication with the 
panel was conducted through an SFWMD web board, which also was available to the public. Publicly 
noticed meeting dates were posted prior to the meetings. Correspondence between SFWMD staff and the 
peer-review panel was conducted through the web board, and all documents, meeting summaries, and task 
deliverables were posted on the web board. Major discussion topics included the model input data sets, 
conceptual model, calibration strategy, model calibration results, and final model documentation. 


2 DESCRIPTION OF THE STUDY AREA 
2.1 Rainfall 


Rainfall is an important part of the hydrologic cycle. The spatial and temporal distribution of rainfall 
influences other variables in the LWCSIM. Long-term (1965 to 2013) historical annual rainfall averages 
54 inches in the LWC Planning Area, and the distribution is shown in Figure 2-1. The driest year recorded 
is 2000, with 41 inches of rainfall, while the wettest year recorded is 1995, with 71 inches of rainfall. 


Rainfall data for the transient simulation period (1999 to 2014) were derived from a combination of Next 
Generation Weather Radar (NEXRAD) and DBHYDRO (the SFWMD’s corporate database) gauge data 
coupled with a simple spatial interpolation technique (Brown 2014). NEXRAD data were used from the 
beginning of May 2002, while gauged data were used from January 1999 through April 2002. NEXRAD 
data provide complete spatial coverage of rainfall amounts using a predetermined grid resolution 
(2 kilometers by 2 kilometers). Four NEXRAD sites operated by the National Weather Service cover the 
LWCSIM domain: KBYX in Key West, KAMX in Miami, KMLB in Melbourne, and KBTW in Tampa. 
Although these sites are not within the LWCSIM boundary, their coverages encompass the model domain. 
The monthly distribution of rainfall across the LWCSIM domain for 1999 and an average of 1999 to 2014 
is shown in Figure 2-24. 
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Figure 2-1. Annual average rainfall distribution in the Lower West Coast Planning Area from 1965 to 
2013. 
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Figure 2-2. Monthly rainfall distribution for 1999 and the average monthly rainfall from 1999 to 2014. 


21.1 Adjustment of NEXRAD Rainfall Data 


A Statistical analysis comparing the daily rain gauge timeseries values with nearby NEXRAD rainfall data 
revealed a seasonal bias in the NEXRAD daily data. The relative bias between the monthly and yearly 
recorded rain gauge and estimated NEXRAD rainfall totals was calculated as follows: 


У NRDi-Y. всі 


Relative bias = - 
УЕСЇ 


x 100 (1) 


Where NRD; and КС; represent the mean rainfall estimated by NEXRAD and recorded by the rain gauge, 
respectively. The mean rainfall difference (RD) can be calculated as follows: 


RD = NRDi — RGi (2) 


The bias was determined to represent overestimation of small rainfall accumulations (< 0.1 inches/day) and 
an underestimation of more extreme events (> 20 inches/day). A method was developed and used to adjust 
the NEXRAD data to reduce this bias. The correction involves identifying a multiplier to apply to each 
daily NEXRAD data value in a manner that a) does not significantly increase or decrease the overall volume 
of water for each year, and b) maintains the same general daily spatial distribution or rainfall patterns. The 
correction does not significantly increase or decrease the annual total NEXRAD data because the NEXRAD 
and rain gauge annual totals have a very good correlation and are very similar. 


To eliminate the bias for other models, the NEXRAD data have been adjusted using a single non-linear 
equation or performing separate adjustments for wet season and dry season with different equations. 
However, the spatial extents of those models were much smaller than this LWCSIM effort required, and 
developing a suitable equation was unnecessary. 


Linear regression of yearly, monthly, and daily rain gauge and NEXRAD rainfall totals was performed to 
evaluate the overall relationship between the variables. Furthermore, the Kolmogorov-Smirnov statistical 
test for normality showed a positive skew of the daily values, indicating a considerable number of small 
values near or equal to zero, which is not surprising because rainfall does not occur every day. Therefore, 
the relationship between NEXRAD data and individual rain gauges was evaluated using the Kendall’s tau 
(t) rank-based correlation method (Helsel and Hirsch 1992). Monthly correlation for each rain gauge and 
NEXRAD pixel pair was calculated using daily values (n > 14), and correlation coefficients were used to 
discern which rain gauges displayed a good relationship (т > 0.75) with the NEXRAD pixels. Rain gauges 
showing a poor relationship with the NEXRAD pixels were excluded from further analysis. Yearly and 
monthly rain gauge and NEXRAD totals were calculated. The bias between the rain gauge and NEXRAD 
data was determined by the ratio between the two (Steiner et al. 1999, Goudenhoofdt and Delobbe 2009, 
Smith and Rodriguez 2017). The resulting value was used as the NEXRAD adjustment factor, calculated 
as follows: 


NRDi 


RGi 
T4 


T4 


AdjF (bias) = X — /Y (3) 
Where AdjF is the multiplicative factor used to adjust the daily NEXRAD value, and RG; and NRD; are the 
monthly rainfall recorded by the rain gauge and estimated by NEXRAD, respectively, during an equally 


match-paired number of days (п > 14). 


The statistical analysis performed did not consistently show a strong correlation when comparing daily rain 
gauge data with NEXRAD data at nearby pixels. However, there was a strong correlation, for most rain 
gauges, when comparing the same data for annual totals and a very good correlation with monthly totals. 


The method developed to reduce bias uses monthly adjustment factors at each rain gauge with a suitable 
correlation and number of observations within the selected month. The adjustment factor is calculated as 
the ratio of the monthly sum of rain gauge values divided by the sum of the nearest monthly NEXRAD rain 
pixel values. Rain gauge stations with adjustment factors were filtered prior to creating NEXRAD bias 
corrections. Included stations must have adjustment factors between 0 and 10, Kendall t correlation values 
less than 0.6, and 15 or more daily values per month. 


To prevent overcorrection of NEXRAD daily values and to include most adjustment factors, the adjustment 
factor was arbitrarily constrained to between 0.7 and 1.3. Adjustments less than 0.7 were assigned the 
minimum adjustment factor of 0.7 and adjustments greater than 1.3 were assigned 1.3. Rain gauges with 
monthly adjustment factors equal to 0 were set to 1 providing no adjustment. Visual examination of the 
resulting bias showed several rain gauges located less than 2 kilometers from each other, which can produce 
inconsistencies in the bias arising from high and low values in spatially close stations. Therefore, the 
resulting monthly bias value at each rain gauge station was normalized using a complete hierarchical cluster 
analysis (Miillner 2013, Maechler et al. 2018), in which bias points are the result of the averaged bias of 
each rain gauge within 2 kilometers. The resulting monthly bias was used as the multiplicative adjustment 
factor to correct the daily NEXRAD pixel value using geostatistical processes. 


Adjustment factors were interpolated across the entire NEXRAD grid using the automated autoKrige 
method described by Hiemstra et al. (2009). The autoKrige function tests different interpolation models 
(linear, spherical, exponential, Gaussian, and Matern [M. Stein’s parameterization]) by estimating 
semi-variograms and selecting the best fit kriging model (Hiemstra et al. 2009). This process produced a 
monthly bias grid with adjustment factors at each NEXRAD pixel. The adjustment factors were applied to 
the uncorrected NEXRAD values to produce corrected NEXRAD values. Point values for the bias and the 
uncorrected and corrected NEXRAD were converted into a raster for visual evaluation. During this process, 
a few inconsistencies associated with the NEXRAD and rain gauge pairing were observed. The 
inconsistencies resulted from mismatched NEXRAD and rain gauge daily values, and the data were 
examined and removed from daily values. Subsequently, the entire process was repeated and a new set of 
raster plots was created. The difference or change in rainfall was calculated by subtracting the corrected 
from the uncorrected NEXRAD, and the results were displayed in a separate raster. The model cell locations 
were used to extract values from the corrected NEXRAD raster and exported as data for each model mesh 
cell. All statistical and spatial data analyses were performed using R programming language (version 3.5.1; 
R Development Core Team 2018). 


The results obtained by applying adjustment factors for each temporal aggregation to make bias adjustments 
and recalculating monthly and annual sums near observed locations were reviewed to ensure the correction 
did not create or eliminate significant amounts of rainfall. Monthly totals provided significantly better 
correction than wet season/dry season totals. There was not a significant benefit to using bi-weekly or 
weekly corrections because the observed and NEXRAD values did not correlate well at those levels of 
discretization. More detailed description of the adjustment of the NEXRAD rainfall data is provided in 
Appendix A. 


The steady-state model used rainfall data from 1999. The maximum monthly rainfall in 1999 (12.9 inches) 
occurred in June, and the minimum amount of rainfall (0.7 inches) occurred in February. Average annual 
rainfall across the model domain in 1999 was 57.3 inches. Spatial distribution of adjusted NEXRAD rainfall 
data for 2005 (wet year) and 2007 (dry year) within the model domain are shown in Figures 2-3 and 2-4, 
respectively. Average rainfall across the model domain was 63.8 and 42.7 inches in 2005 and 2007, 
respectively. 
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Figure 2-3. Rainfall distribution across the model domain in 2005 (wet year). 
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Figure 2-4. Rainfall distribution across the model domain in 2007 (dry year). 
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2.2 X Evapotranspiration 


Evapotranspiration (ET) is the sum of the evaporation from water bodies and transpiration losses from plant 
systems to the atmosphere. ET plays a major role in the hydrologic cycle in the LWC Planning Area because 
the water table is near land surface throughout the year. In fact, ET causes the largest loss of water from the 
system. 


Reference ET is the rate of ET from a hypothetical reference crop with an assumed crop height of 
0.12 meters, a fixed surface resistance of 70 seconds per meter, and an albedo of 0.23, closely resembling 
the ET from an extensive surface of green grass of uniform height, actively growing, well-watered, and 
completely shading the ground (Jensen et al. 1990, Doherty et al. 2010). Potential ET is the amount of water 
used by a crop and can be estimated by multiplying the reference ET by a crop water use coefficient (K.). 
The K. is specific to a particular crop or vegetation and can vary from month to month. 


Reference ET data were obtained from the SFWMD’s updated data set (Brown 2013). The updated data 
have a domain consistent with the jurisdictional boundary for the SFWMD and were created from two data 
sets. The first data set was the 51-year hydrological re-analysis (HYDROS1) and the second data set was 
generated by the North American Land Data Assimilation System (NLDAS). The long-term data set was 
created incrementally, and the completed set is believed to be the best available because the data from 1979 
to present used the NLDAS-based satellite data (Brown 2013). 


Historical reference ET distribution from 1965 to 2014 is shown in Figure 2-5. The average value over this 
period is 57 inches/year. The monthly distribution of reference ET across the LWCSIM domain is shown 
in Figure 2-6. The maximum reference ET rate (6.24 inches) occurred in April, and the lowest reference 
ET rate (3.01 inches) occurred in December (Smajstrla 1990). Figure 2-6 also includes the average 
reference ET rates for the LWCSIM transient calibration period (1999 to 2014). Figures 2-7 and 2-8 show 
the spatial distribution of reference ET across the LWCSIM domain for 2005 (wet year) and 2007 (dry 


year). 


11 


70 
60 
50 
_ 40 
oO 
Ф 
> 
У 
[^4] 
o 30 
с 
o 
Е 
20 
10 
0 
o A 9 о AV © ау o> do «А OD дуо? о SN 0 су > A OQ лу к 
РЯ 
Figure 2-5. Annual average reference evapotranspiration from 1965 to 2013. 
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Figure 2-6. Monthly reference evapotranspiration distribution for 1999 and average monthly reference 


evapotranspiration distribution from 1999 to 2014. 
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Figure 2-7. Spatial distribution of reference evapotranspiration in 2005 (wet year). 
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Figure 2-8. Spatial distribution of reference evapotranspiration in 2007 (dry year). 
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2.3 Весһагде 


Within Ше LWC Planning Area, groundwater recharge primarily occurs through rainfall and, to a lesser 
extent, through irrigation return flows. Groundwater flow occurs from high to low elevation areas, the 
amounts of which can vary dramatically throughout the year. 


Although the average annual rainfall and temperatures are relatively similar across the planning region, 
large variations in spatial and temporal distributions occur daily. These differences are obtained daily from 
rain gauges, NEXRAD, and climatic stations. This site-specific information is used when calculating the 
supplemental irrigation demands using the Agricultural Field-Scale Irrigation Requirement Simulation 
(AFSIRS) program (Smajstrla 1990) for specific properties. 


The amount of rainfall available for recharge to the SAS is reduced by runoff and ET. In this area, ET 
generally accounts for 60% to 70% of rainfall, and the remaining 30% to 40% runs off the land into drainage 
networks/streams/rivers or percolates into the ground as recharge. 


The amount of water available from the SAS as recharge to the IAS depends on several variables. Areas 
where the soil is well drained are potentially good areas of moderate to high recharge capacity, while in 
areas with poorly drained soil, the recharge potential could be near zero because of higher runoff and ET in 
these areas. A second variable is the thickness and characteristics of the confining units. Areas in which 
a) the confining units are thin, b) the aquifers have a higher hydraulic conductivity (K), or c) local aquifers 
make up a dominant percentage of the overall unit thickness are prime areas for recharge because flow 
through these units is less restrictive. The third primary variable is the downward hydraulic gradient 
between the SAS and IAS. For recharge to occur, water levels in the SAS need to be higher than water 
levels in the underlying IAS. 


Figures 2-9 and 2-10 show the spatial distribution of annual average gross recharge provided to 
MODFLOW in 2005 (wet year) and 2007 (dry year), respectively. Gross recharge is defined as follows: 


Gross Recharge = Rain + Irrigation — Initial Abstraction - Runoff — Unsaturated/Vadose Zone ET (4) 


Gross recharge was estimated using a pre-processor, which is described in Section 5.7 and should not be 
confused with net recharge, which is the combined effect of applied recharge and ET. The 
ET-Recharge-Runoff pre-processor separates runoff from rainfall using the United States Department of 
Agriculture (1986) National Resources Conservation Service (NRCS) Curve Number (CN) approach and 
estimates the vadose/unsaturated zone ET. Gross recharge is provided as input to MODFLOW. Net 
recharge is defined as follows: 


Net Recharge = Gross Recharge — Saturated Zone/Groundwater ET (from MODFLOW) (5) 


Net recharge is the SAS recharge that remains after MODFLOW rejects some of it as groundwater ET using 
the ET package. Thus, net recharge is calculated after the MODFLOW simulation is complete. The highest 
gross recharge occurs in the southern portion of the model domain, which is primarily wetland areas such 
as Everglades National Park and Big Cypress National Preserve (Figures 2-9 and 2-10). This is because, 
in wetland areas, runoff is routed through the Wetlands package and is not included in the gross recharge 
equation (4). Also, the unsaturated zone is not present in wetland areas during most of the year (shallow 
water table condition); thus, very small unsaturated zone ET is removed from the system. The lowest gross 
recharge is applied to the northeastern portion of the model domain (Hendry and Glades counties), where 
the unsaturated zone is present (deeper water table condition); thus, a significant amount of unsaturated 
zone ET is removed from the system. 
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Figure 2-9. Spatial distribution of gross groundwater recharge rates in 2005 (wet year). 
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Figure 2-10. Spatial distribution of gross groundwater recharge rates in 2007 (dry year). 
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2.4 бой Infiltration Groups 


Infiltration rates of soils in the LWCSIM domain vary widely, from high infiltration capacity (low runoff 
potential) to low infiltration capacity (high runoff potential). The runoff-rainfall potential depends on the 
hydrologic soil group, as classified by the NRCS. Soils are classified into four hydrologic soil groups with 
specific letter designations: A, B, C, and D designations indicate high, moderate, low, and poor infiltration 
capacities, respectively. Dual designations, such as A/D, B/D, or C/D are assigned where a high water table 
may limit infiltration. In these cases, the first letter applies to the drained condition, while the second applies 
to the undrained condition. The ET-Recharge-Runoff pre-processor that generates recharge to the LWCSIM 
uses the hydrologic soil group information during the rainfall-runoff separation process of the NRCS CN 
approach. 


With 59% of areas classified as A/D, soils in the LWCSIM domain predominantly have high infiltration 
rates that may be limited by a shallow water table condition (Figure 2-11). Most A/D soils in southern 
Collier County behave as D soils due to the undrained condition (shallow water table). The next most 
prevalent soil group is C/D (16%); these are poorly drained, low-infiltration soils commonly observed in 
the central and northwestern areas of the LWCSIM domain. The third most common soil group is B/D 
(11%), which consists of moderately drained, somewhat low-infiltration capacity soils found throughout 
the model domain. 
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Figure 2-11. Spatial distribution of hydrologic soil groups across the model domain, as classified by the 
Natural Resource Conservation Service. 
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2.5 Land Use 


Land use represents a spatial distribution of pervious and impervious surfaces used to separate runoff and 
infiltration from total rainfall and irrigation. The steady-state model used the SFWMD's 2000 land use map, 
which was produced using data obtained between 1998 and 2000. The transient version of the model used 
four land use maps: 1) 2000, representing 1999 to 2001; 2) 2004, representing 2002 to 2006; 3) 2010, 
representing 2007 to 2011; and 4) 2014, representing 2012 to 2014. These maps used the Florida Land Use 
and Cover Classification System (FLUCCS) (SFWMD 2009), which updates land use every 5 years based 
on the SFWMD’s enterprise data collection. 


Land cover and land use areas primarily were used in the model to estimate maximum ET and recharge 
rates, AG and L/R irrigation water demands, and irrigation return flow rates. The NRCS CN method (United 
States Department of Agriculture 1986) was used in an external pre-processing routine to separate runoff 
from the combination of rainfall and irrigation as a function of land cover and properties of the surface soil. 
Land use-based crop type and permeability were developed for the FLUCCS codes. 


Table 2-1 provides the various land use classifications and acres within the LWCSIM domain. 
Figures 2-12 through 2-15 depict the land use/land cover maps for 2000, 2004, 2010, and 2014, 
respectively. 


Table 2-1. Land use summary for 2000, 2004, 2010, and 2014 in the Lower West Coast Planning Area. 


Land Use 2000 (acres) 2004 (acres) 2010 (acres) 2014 (acres) 

Barren Land 24,919 26,676 29,336 42,327 
Citrus 230,519 223,901 197,053 173,996 
Forest 298,340 284,009 309,884 297,247 
Golf 17,511 23,097 24,533 25,347 
Improved Pasture 391,506 370,311 397,770 389,238 
Non-Urban Natural Land 311,619 266,190 227,613 261,059 
Nursery 4,010 4,312 8,787 8,395 
Pasture 123,174 109,669 94,552 101,046 
Row Crops 102,465 103,384 107,549 79,508 
Sod 1,262 684 1,179 2,326 
Sugar 183,730 190,867 187,823 203,918 
Urban 232,584 308,853 374,024 372,426 
Water 501,508 458,332 485,624 447,284 
Wetlands 1,721,845 1,724,757 1,694,417 1,697,201 
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Figure 2-12. 2000 land use/land cover map. 
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Figure 2-13. 2004 land use/land cover map. 
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Figure 2-14. 2010 land use/land cover map. 
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Figure 2-15. 2014 land use/land cover map. 
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3 PHYSICAL SETTING 


3.1 Surface Water Flow System 


Primary surface water sources in the LWCSIM domain include the Caloosahatchee River (C-43 Canal) and 
connected water bodies, such as the Townsend Canal, Roberts Canal, and City Ditch. Cape Coral and Big 
Cypress Basin (BCB) canal systems also provide surface water supply, and to a lesser extent, local storm 
water ponds provide water for landscape irrigation. AG is the primary consumer of surface water in the 
LWC Planning Area. 


The BCB canal system in Collier County was constructed as a surface water drainage system. Since 2000, 
improvements to structures, operations, management, and monitoring have resulted in an estimated 
850 acre-feet of additional surface water storage in BCB canals (SFWMD 2012). The BCB capital 
improvement program includes projects for the Golden Gate canal system, Henderson Creek, FAKA Union, 
and Barron River. These projects and others provide water resource benefits through reduction of 
over-drainage and restoration of groundwater and surface water levels to more natural conditions. In 
addition to providing environmental benefits, these improvements enhance water supply opportunities by 
increasing groundwater storage and improving the timing and duration of surface water discharges. As a 
result, the canal system holds more water for longer periods of time than previously capable, capturing 
water that used to be lost to tide. 


The City of Naples Utility Department is developing surface water sources from the BCB canal system to 
supplement its reclaimed water supply. Cape Coral Utilities also uses a freshwater canal system to augment 
the City of Cape Coral’s reclaimed water supply for residential and commercial landscape irrigation. The 
City of Marco Island (Marco Island Utilities) uses surface water from Henderson Creek/Marco Lakes, and 
Lee County Utilities uses some surface water from the Caloosahatchee River. 


32 Hydrostratigraphy 


3.2.1 Geological Framework 


South Florida, including the LWC Planning Area, is underlain by a thick sequence of carbonate and clastic 
sedimentary rocks ranging in age from Paleocene to Recent. Because this study is limited to the SAS and 
IAS, only the deposits ranging in age from late Oligocene through the Holocene are discussed. The 
stratigraphic units of interest to this study include the Peace River and Arcadia Formations of the Hawthorn 
Group, the Tamiami Formation, and the undifferentiated sediments of Holocene/Pleistocene age. 
Figure 3-1 shows the general geology and hydrostratigraphy within the LWCSIM domain. 


25 


System Hydrogeologic Unit Lithostratigraphic Unit 


Undifferentiated Holocene/Pleistocene 
WATER TABLE AQUIFER 


Pinecrest Sand Member 


Bonita Springs Marl Member / 


REE TENS UNIT Caloosahatchee Clay Member 


Tamiami 
Formation 


LOWER TAMIAMI AQUIFER 
UPPER HAWTHORN CONFINING UNIT 

| CLASTIC ZONE 

| CARBONATE ZONE 
MID-HAWTHORN CONFINING UNIT 


Ochopee Limestone Member 


SANDSTONE AQUIFER (SA) Peace River Formation 


Apnis siy} уо 12efqns 


MID-HAWTHORN AQUIFER 
LOWER HAWTHORN CONFINING UNIT 


LOWER HAWTHORN 
PRODUCING ZONE 


Intermediate Aquifer | Surficial Aquifer 


Hawthorn Group 


Arcadia Formation 


UPPER FLORIDAN AQUIFER 


н отет Avon Park Formati 
qe PERMEABLE ZONE — ám 


Floridan Aquifer System 


LOWER FLORIDAN AQUIFER Oldsmar Formation 
Cedar Keys Formation 
| [| SUB-FLORIDAN CONFINING UNIT 


Figure 3-1. Generalized geologic and hydrogeologic units of the Lower West Coast Planning Area 
(From: Geddes et al. 2015). 


The lithology of the undifferentiated surficial Holocene/Pleistocene sediments is highly variable. 
Medium- to fine-grained quartz sand, fossils, clays, and some freshwater limestone and marl are present 
within this unit. These extend to the top of the Tamiami Formation. 


In a few areas, the Tamiami Formation is entirely absent, and the surficial sediments rest directly on top of 
the Peace River Formation. The undifferentiated surficial sediments grade into the Anastasia Formation to 
the east and into the Miami Limestone to the south (Bryan et al. 2013). The Tamiami Formation is composed 
of two units and nine members, none of which are present throughout the entire project area. The upper 
portion of the Tamiami Formation is predominantly marl, and the lower portion is the Ochopee Limestone. 
Lithology of this unit varies from fine- to coarse-grained sand and fossiliferous limestone (Scott 2001). The 
presence of these two units varies spatially, and the Ochopee Limestone is absent in much of southwestern 
Hendry and northeastern Monroe counties. 


The Peace River Formation of the Hawthorn Group is Miocene in age and consists of clays and carbonates 
interbedded with quartz sands. Phosphate may be gravel to sand sized. Approximately two-thirds of the 
formation is siliciclastic, and one-third (typically the lower portion) is carbonate. The Peace River 
Formation underlies the entire project area. 


The underlying Arcadia Formation of the Hawthorn Group is predominantly carbonate and spans the entire 
project area. The contact between the two formations may be distinct or gradational. The Arcadia Formation 
is primarily dolostone and limestone with beds of clay, quartz sand, and phosphate grains (Scott 1988). The 
base of the Arcadia Formation is confining in nature and primarily clay and mud. 
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3.2.2  Hydrostratigraphic Framework 


Due to the presence of a thick clay layer underneath the SAS/IAS systems in the LWC Planning Area, the 
Floridan aquifer system typically is hydraulically disconnected from the overlying SAS/IAS system; 
therefore, the Floridan aquifer system is not included in the LWCSIM and will not be discussed further. 


The hydrogeology of southwestern Florida is complex. Lateral facies changes and variable thicknesses lead 
to large local variations in hydrogeologic units. The heterogeneous nature of the units and the sparse 
availability of data in some places pose difficulties for regional mapping. The hydrogeologic units of 
interest for the LWCSIM effort include the WTA, Tamiami confining unit (TCU), LTA, SSA (clastic and 
carbonate zones), Mid-Hawthorn confining zone, and the MHA. Geddes et al. (2015) provides further 
details on these hydrogeologic units. 


The WTA is composed primarily of quartz sand and shell, with minor amounts of organic material. A dense 
limestone cap rock is present in some areas. The WTA is absent or insignificant in places within the 
LWCSIM domain (e.g., south-central Collier County), and the basal confinement is geographically 
variable. Where present, the Bonita Springs Marl and low-permeability portions of the Pinecrest Sand 
facilitate confinement at the base of the WTA. Confinement between the WTA and the underlying LTA, 
however, is not consistent. Where the TCU is absent or insignificant, the LTA, by definition, could be 
considered part of the WTA. 


The LTA is predominantly a sandy, biogenic limestone and calcareous sandstone. It encompasses all the 
water-producing limestone and, in some areas, portions of the underlying permeable sand. The upper 
confinement (TCU) is absent or insignificant in several areas throughout the LWCSIM domain. In northern 
portions of the modeled area (i.e., Charlotte County into the Southwest Florida Water Management 
District), well reports typically do not distinguish sub-units of the SAS. Throughout most of the LWCSIM 
domain, the less permeable clay and fine-grained sands of the Peace River Formation provide basal 
confinement (Upper Hawthorn confining unit) to the LTA. However, in some areas, this confinement is 
absent or insignificant. In some locations, the LTA or undifferentiated SAS lay directly on top of the SSA. 


The SSA, except where absent or insignificant, is contained entirely within the Peace River Formation of 
the Hawthorn Group and is part of the IAS. It typically occurs as two distinct permeable units: an upper 
clastic zone and a lower carbonate zone. The SSA is composed of sandstones, sandy limestones, dolostones, 
and calcareous sands. These may be contiguous or separated by varying amounts of low-permeability silt 
and clay. Where a confining unit is present, the SSA is separated from the LTA by lower-permeability clays 
and dolosilts of the Peace River Formation (Upper Hawthorn confining unit). The SSA is separated from 
the underlying MHA by low-permeability clays and marls of the basal Peace River Formation 
(Mid-Hawthorn confining unit), which are present throughout the LWCSIM domain. The SSA is absent or 
insignificant throughout the eastern portions of the model domain and some smaller areas along the Gulf 
coast. 


The MHA, except where absent or insignificant, is composed of biomicritic limestone, phosphate, shell, 
and lime mud. It lies entirely within the Arcadia Formation of the Hawthorn Group. The MHA is separated 
from the overlying SSA by the low-permeability clays and marls of the basal Peace River Formation 
(Mid-Hawthorn confining unit). Where the SSA is absent or insignificant, the entire thickness of the Peace 
River Formation isolates the MHA from the overlying SAS. The confinement from the underlying Lower 
Hawthorn producing zone (not simulated in this model) consists of carbonate muds and terrigenous clays 
of the Upper Arcadia Formation (Lower Hawthorn confining unit) and is present throughout the study area. 
For the most part, use of the MHA for water supply purposes occurs in the western part of the LWCSIM 
domain. Some investigators have defined the MHA as two distinct zones: a shallow zone occurring 
primarily in the Fort Myers-Cape Coral area, and a deeper zone occurring primarily in the Estero-Bonita 
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Springs-Fort Myers Beach region. Lacking the data to support а more detailed breakdown of the aquifer 
throughout the LWCSIM domain, the MHA is considered a single unit for this modeling effort. 


3.2.8 Conceptualized Groundwater Flow 


Regional groundwater flow is consistent across the WTA, LTA, and SSA. Loosely following land surface 
elevations, regional flow patterns in the LWCSIM domain can be subdivided into two components. In the 
northern portion of the modeled area, flow is in a general northeast-to-southwest direction from a regional 
high along the southern portions of the Highlands Ridge in Polk and Highlands counties towards the Gulf 
coast and the Caloosahatchee River. In the central portion of the modeled area, flow in all three aquifers 
generally radiates out in all directions from a regional high located close to the Hendry-Collier county 
border north of the Town of Immokalee. 


Regional flow patterns in the MHA show a general northeast-to-southwest flow pattern throughout the 
LWCSIM domain. Due to its greater degree of confinement and lower levels of productivity, the MHA is 
more likely to show impacts from well pumping than the shallower, more productive aquifers. An area of 
regional decline in the potentiometric surface of the MHA is centered around Cape Coral, which 
corresponds to an area of historically high pumpage from PS wells and small domestic and irrigation wells. 


3.2.4 Data Collection and Analysis 


Data collection and assembly of the corresponding computer-model input data sets are among the most 
time-consuming elements of model development, particularly for a model with as large a geographic extent 
and complex hydrostratigraphy as the LWCSIM. The following sections describe the data used to construct 
the steady-state model and how the data were obtained. 


3.2.5  Hydrostratigraphic Data 


The conceptual groundwater model was divided into hydrostratigraphic units following the primary 
water-producing units and confining units identified by Geddes et al. (2015). Data from their work and the 
Districts DBHYDRO database provided most of the information used to construct the conceptual model. 
The groundwater flow model was constructed based on hydrostratigraphic units and may not necessarily 
follow geological formation contacts. 


Once all data were assimilated and checked for quality assurances, they were combined into a single data 
set to develop the model layers. To solve the flow equations in MODFLOW 2005, the continuous presence 
of aquifers is required. In other words, aquifer or layer thickness must be greater than zero everywhere in 
the model domain regardless of whether or not the aquifer is present. To satisfy this condition and to ensure 
numerical stability during model execution, a minimum thickness of 30 ft was assigned to all layers. Where 
aquifers are absent, the hydraulic properties of the aquifer(s) below were used. The hydrostratigraphic 
layers were created by kriging the surfaces of each major hydrostratigraphic unit. Appendix B, Figures B-1 
to B-10 depict the surface elevations for each of the nine model layers, beginning with the top of layer 1 
(the topographic surface), then progressing downward through each layer's bottom surface. The bottom of 
a layer represents the top of the subjacent layer. 


Some of the hydrostratigraphic control points indicate a unit may be extremely thin or absent (defined as 
5 ft or less). Appendix B, Figures B-11 to B-19 are isopach maps of the aquifers and confining zones. 
Simulated layer top elevations are show in Appendix B, Figures B-20 to B-28. Appendix B, Figures B-29 
to B-37 provide the simulated layer thicknesses. 
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3.2.6 Hydraulic Data 


Multiple types of hydraulic data are required to develop a multi-layered groundwater flow model. For this 
section of the model development discussion, the data are divided into the following four categories: 


Hydraulic conductivities and storage coefficients 
Pumpage rates 

Boundary conditions 

Initial conditions 


The first two categories are directly related to aquifer properties and hydraulic data, while the remaining 
two represent aquifer stresses and observed conditions. 


Data collection procedures generally require that data be collected both from within and outside the model 
domain to define conditions along the model boundaries as accurately as possible. The discussion that 
follows focuses on the data collected within the model boundaries unless specifically stated otherwise. 


3.2.7 Hydraulic Conductivities/Storage Coefficients 


Hydraulic conductivity (K) is an important parameter in development of a groundwater flow model. It 
represents the aquifer’s ability to transmit water under hydraulic gradients and hydraulic stresses such as 
pumping wells. When multiplied by the aquifer thickness, the resulting term is called transmissivity, which 
can be readily obtained from aquifer performance tests (APTs). For the LWCSIM, aquifer and confining 
unit top and bottom elevations are considered static input parameters. The thickness of the aquifer, and 
consequently the transmissivity of the aquifer, was calculated internally by the model code. Therefore, 
vertical and horizontal hydraulic conductivities were used in place of transmissivity values. Figures 3-2 to 
3-10 present the interpolated horizontal hydraulic conductivity (Kh) surfaces for each model layer used to 
define this value in each model cell. The figures were created by interpolating the actual field data values 
collected from APTs, specific capacity, and slug tests and then using kriging and spline with barriers. When 
multiple test data points were available in close proximity, priority was given to APT values due to their 
inherent higher quality. The surfaces were used as the starting point for the steady-state model calibration. 
The K distribution in all nine layers was adjusted as a result of steady-state and transient calibration. The 
adjusted K distributions are presented and discussed in the Calibration section. A discussion of how layer 
pinch-outs (absences) were represented within the K distribution of each layer also is presented in 
Section 6. 


Another important parameter needed for the transient simulations is the storage coefficient of the aquifers. 
Storage coefficients can be readily obtained from APTs. Like the relationship between K and transmissivity, 
the model code requires a specific storage value, which is the quotient of the aquifer storativity divided by 
aquifer thickness. Storage coefficient (storativity) values were entered directly into each model layer as a 
property based on APT data. 
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Figure 3-2. Horizontal hydraulic conductivity of the Water Table aquifer. 
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Figure 3-3. Horizontal hydraulic conductivity of the Tamiami confining unit. 


31 


SWFWMD 


1 Esoto 


LAKE 
OKEECHOBEE 


GULF OF 
MEXICO 


Kx from Field Data 


Layer 3 (Lower Tamiami 
Aquifer) 


Test 
© APT ft/day 
© SC ftday 
Aquifer Absence 
Feet 
ШИШ о-5 
Hydraulic Conductivity 
(#аау) 
ШШ 0.01 - 0.1 
EB 01-1 


E 1-10 KOT [J LWCSIM Model Extent 
Г] 10 - 100 SFWMI LWCSIM Active Model Boundary 


| 100 - 1000 { 15 20 Мйез 
ШЕ 1000 - 10000 Е 30 Kilometers 


\\ad.sfwmd gov \\ Е Figures\mxd\KxfromF ieldData.mxd 


Figure 3-4. Horizontal hydraulic conductivity of the Lower Tamiami aquifer. 
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Figure 3-5. Horizontal hydraulic conductivity of the Upper Hawthorn confining unit. 
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Figure 3-6. Horizontal hydraulic conductivity of the Sandstone aquifer — clastic zone. 


34 


SWFWMD 


SARASOTA 


LAKE 
OKEECHOBEE 


GULF OF (ҚР АЛАЛА 


MEXICO 


Kx from Field Data 


Layer 6 (Sandstone Confining 
Unit) 


Test 
O APT ft/day 
© SC Жадау 
Aquifer Absence 
Feet 
ШИШ о-5 
Hydraulic Conductivity 
(#аау) 
ШШ 001-01 y - 
ШЕН 0.1-1 | p 


EH 1-10 У | E LWCSIM Model Extent 
Г 1] 10 - 100 ИМ! ГД LWCSIM Active Model Boundary 


[SH 100 - 1000 — | x в юм 
ШШ 1000 - 10000 2 | Тукай 


“MIAMI-DADE - 


- - — m RU. 


BAS 


\\ad.sfwmd.gov \dfsroot\GS\GSBiz\WS\LWCSIM\Report_Figures\mxd\KxfromF ieldData.mxd 


Figure 3-7. Horizontal hydraulic conductivity of the Sandstone aquifer confining unit. 
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Figure 3-8. Horizontal hydraulic conductivity of the Sandstone aquifer — carbonate zone. 
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Figure 3-9. Horizontal hydraulic conductivity of the Mid-Hawthorn confining unit. 
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Figure 3-10. Horizontal hydraulic conductivity of the Mid-Hawthorn aquifer. 
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3.3 Water Use 


Creation of the MODFLOW WEL package for the LWCSIM required extensive data related to historical 
water use. The overall goal was to use the best available information related to historical use as well as the 
best available methods for estimating historical use when no withdrawals were reported. Therefore, the data 
used for the WEL package were extracted from individual water use permits, actual water use reported to 
the SFWMD, and demand estimates. 


3.3.1 Water Use Summary 


There are six water use categories used in water supply planning in South Florida: Public Supply (PS), 
Agriculture (AG), Landscape/Recreational (L/R), Commercial/Industrial/Institutional (СП), Domestic 
Self-Supply (DSS), and Power Generation (PG). There are no Power Generation permits within the 
SAS/IAS, so the category is not discussed in this report. 


1999 Steady-State Water Use Summary 


Table 3-1 shows the total SAS/IAS water use demands in 1999 for each county and water use category 
within the LWC Planning Area. Figure 3-11 shows the total amount and corresponding percentage for each 
water use classification. AG is the largest water use category within the model domain, making up 71% of 
total water use. Figure 3-12 shows the total amount and corresponding percentage for each county within 
the model domain. Collier County is using the largest volume of water. Figure 3-13 shows the distribution 
of water withdrawals from each layer of the model. The LTA is the most used aquifer in the LWCSIM 
domain, with approximately 228 mgd withdrawn. 


Table 3-1. Water use, in million gallons per day, by county and use type, of the surficial and 
intermediate aquifer systems, for the Lower West Coast Planning Area in 1999. 


County AG CII DSS L/R PS 
Charlotte 9.12 0.00 0.03 0.00 0.05 
Collier 126.95 1.13 2.26 28.71 63.85 
Glades 9.46 0.00 0.00 0.10 0.01 
Hendry 177.30 0.71 0.00 0.51 0.00 
Lee 63.93 0.25 15.13 36.09 11.09 


AG = Agriculture; СП = Commercial/Industrial/Institutional; DSS = Domestic Self-Supply; L/R = Landscape/Recreational; 


PS = Public Supply. 
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Figure 3-11. Water use distribution by use category, in million gallons per day, of the surficial and 
intermediate aquifer systems, for the Lower West Coast Planning Area in 1999. 
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Figure 3-12. Water use distribution by county, in million gallons per day, of the surficial and 
intermediate aquifer systems, for the Lower West Coast Planning Area in 1999. 
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Figure 3-13. Water use distribution by aquifer layer, in million gallons per day, of the surficial and 
intermediate aquifer systems, for the Lower West Coast Planning Area in 1999. 


Transient Period (1999 to 2014) Water Use Summary 


Table 3-2 shows the water use demands, by use type, for the transient calibration and verification periods 
(1999 to 2014). Figure 3-14 shows the trends in water use demands from 1999 to 2014. AG demands were 
consistently higher than the other water use categories throughout the calibration and verification periods. 
The graph also shows that L/R demands increased in 2007. DSS demand increased gradually up through 
approximately 2007. 


Table 3-2. Water use demands by use type from 1999 to 2014, in million gallons per day, of the surficial 
and intermediate aquifer systems, for the Lower West Coast Planning Area. 


Year AG СП DSS L/R PS 

1999 387.49 2.08 17.42 68.17 11.15 
2000 447.99 1.97 19.88 70.81 11.43 
2001 389.13 1.66 23.29 67.33 34.64 
2002 392.03 2.19 25.77 63.79 38.54 
2003 318.63 1.18 30.41 92.56 36.17 
2004 378.47 1.15 36.91 104.47 43.40 
2005 332.45 1.24 43.40 93.72 64.73 
2006 426.60 1.05 51.86 112.27 71.30 
2007 428.08 0.96 55.33 112.09 75.74 
2008 488.67 2.22 54.79 130.99 69.84 
2009 538.18 3.42 55.11 136.12 69.80 
2010 433.89 4.28 56.78 120.87 64.73 
2011 514.23 4.88 51.51 130.89 64.49 
2012 512.72 4,75 51.72 126.38 63.43 
2013 447.07 4.06 51.88 125.78 61.22 
2014 425.25 3.53 52.26 118.06 62.04 


AG = Agriculture; СП = Commercial/Industrial/Institutional; DSS = Domestic Self-Supply; L/R = Landscape/Recreational; 
PS = Public Supply. 
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Figure 3-14. Water use summary for the Lower West Coast Planning Area, 1999 to 2014. 


The following sections describe how the various water use information was compiled and incorporated into 
the model. 


3.3.2 Public Supply 


PS is the second largest water use category within the LWCSIM domain, with 32 individual water use 
permits and 509 facilities. These permits and facilities are spread throughout Charlotte, Glades, Lee, Collier, 
and Hendry counties, with most facilities located in Lee and Collier counties. Figure 3-15 shows the 
location of PS pumping facilities within the LWCSIM domain. 


Historical withdrawal records for PS are available throughout the LWCSIM domain and the calibration 
period via reported values submitted by permittees to the SFWMD's Water Use Bureau. Raw water 
volumes pumped during a given month and year are reported by the corresponding utility. Throughout the 
calibration period, individual well pumping records were reviewed and used when available. Historical 
individual well pumping records were used to create facility distribution ratios. Facility distribution ratios 
were developed based on per-well historical usage and applied to the period of time when per-well pumpage 
records were not available. At the beginning of the simulation period, historical pumping sometimes was 
only reported by permit (a summation of all wells providing water for the utility) or by aquifer (a summation 
of all wells withdrawing water from a specific aquifer). In these cases, it was assumed that facility 
distribution ratios from later years were applicable in preceding years. 


Production wells with open-hole sections penetrating multiple layers were identified and examined. For 


these wells, the withdrawal at the affected well location was divided among the productive layers (excluding 
confining units) proportional to the transmissivities of the contributing open intervals. 
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Figure 3-15. Public Supply pumping facilities using the surficial and intermediate aquifer systems within 
the Lower West Coast Planning Area. 
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3.3.3 Agriculture and Landscape/Recreational 


AG and L/R combined are the largest water use in the LWC Planning Area. These water uses are primarily 
in north-central Collier County, eastern Lee County, Hendry County, and Glades County. The major 
agricultural crops are citrus, sugar cane, turf, improved pasture, and various vegetable crops (e.g., beans, 
tomatoes, melons, other row crops). AG can be subdivided into agriculture, diversion and impoundment, 
aquaculture, livestock, and nursery, while L/R can be subdivided into golf and landscape. The crop type for 
L/R is turf. 


Within the LWCSIM domain, there are 4,504 water use permits and 11,563 facilities in the AG and L/R 
categories. These permits and facilities are spread throughout Lee, Collier, Hendry, Glades, and Charlotte 
counties, with most permittees and facilities located in Lee, Collier, and Hendry counties. Figure 3-16 
shows the locations of AG and L/R pumping facilities within the LWCSIM domain. 


Throughout the calibration period, pumping records by permit and by individual well were reviewed. 
Unfortunately, withdrawal records for the calibration period for AG and L/R users was very limited in 
availability. Therefore, irrigation demands were determined using the AFSIRS program (Smajstrla 1990). 
AFSIRS provides an estimate of daily irrigation requirements based on observed rainfall and ET rates, crop 
types, and land use. A detailed description of how AFSIRS is used to estimate AG demand is described at 
the end of this section. The SFWMD’s water use permit database was queried to determine the predominant 
crop type and irrigation efficiency for each permit. This information was input into AFSIRS, and irrigation 
requirements were calculated for each day of the simulation period for each individual water use permittee. 
Demands were summarized to produce a yearly demand in million gallons per year. Results from AFSIRS 
for AG and L/R permittees are summarized in Tables 3-3 and 3-5, respectively. Tables 3-4 and 3-6 provide 
the acreages by county for AG and L/R uses, respectively. Table 3-7 shows the acreage by crop type for 
each of the land use/land cover maps. During the years with significant increases in acreage, there also are 
corresponding increases in demand. The biggest AG users in the LWCSIM domain are Babcock Ranch, 
located in Charlotte County, and Southern Division Ranch Unit 1, located in Hendry County. 
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Figure 3-16. Agriculture and Landscape/Recreational pumping facilities using the surficial and 
intermediate aquifer systems in the Lower West Coast Planning Area. 
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Table 3-3. Agriculture water use demand by county, in million gallons per day, of the surficial апа 


intermediate aquifer systems within the Lower West Coast Planning Area. 


Year Charlotte Glades Lee Hendry Collier 
1999 23.52 4.09 29.32 135.30 115.86 
2000 34.18 7.97 37.57 194.14 142.98 
2001 30.76 6.44 38.72 188.89 124.85 
2002 26.50 6.76 28.37 140.71 111.34 
2003 30.54 6.37 26.94 140.62 105.51 
2004 35.45 7.89 31.43 163.38 129.86 
2005 33.17 7.02 30.01 146.83 121.06 
2006 37.90 9.18 35.91 192.45 143.54 
2007 29.19 10.38 31.44 171.63 136.63 
2008 23.76 7.56 24.61 150.75 114.63 
2009 30.79 11.18 33.94 202.72 153.81 
2010 26.96 8.58 24.72 151.28 105.11 
2011 31.05 9.86 27.90 184.44 139.14 
2012 27.29 8.63 26.77 142.77 134.02 
2013 24.22 7.90 24.66 144.27 122.37 
2014 22.98 6.70 21.14 126.85 108.59 


Table 3-4. Agriculture acreage by county within the Lower West Coast Planning Area and model 


domain. 
Year Charlotte Glades Lee Hendry Collier 
2000 15,755 42,799 25,841 240,669 68,020 
2004 15,493 43,184 23,931 229,687 68,683 
2010 11,081 47,347 20,418 237,883 65,475 
2014 11,346 45,464 17,068 216,104 63,725 


Table 3-5. Landscape/Recreational water use demand by county, in million gallons per day, of the 
surficial and intermediate aquifer systems within the Lower West Coast Planning Area. 


Year Charlotte Glades Lee Hendry Collier 
1999 0.00 0.01 19.43 0.25 17.36 
2000 0.00 0.01 21.77 0.29 18.17 
2001 0.00 0.01 20.60 0.26 15.87 
2002 0.00 0.01 28.28 0.14 20.90 
2003 0.00 0.01 27.21 0.14 19.13 
2004 0.00 0.01 28.34 0.17 21.87 
2005 0.00 0.01 21.94 0.14 18.20 
2006 0.00 0.01 31.06 0.20 26.81 
2007 0.00 0.13 45.88 0.99 35.13 
2008 0.00 0.09 39.68 0.79 29.61 
2009 0.00 0.13 46.06 1.05 33.80 
2010 0.00 0.10 37.10 0.87 24.92 
2011 0.00 0.12 43.49 1.01 31.37 
2012 0.00 0.12 48.42 1.05 32.16 
2013 0.00 0.12 45.20 1.00 33.64 
2014 0.00 0.10 40.49 0.89 29.88 
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Table 3-6. Landscape/Recreational acreage by county within the Lower West Coast Planning Area. 


Year Charlotte Glades Lee Hendry Collier 
2000 21 477 13,059 1,413 13,041 
2004 21 422 19,750 908 17,593 
2010 25 552 25,432 2,007 21,259 
2014 76 357 28,775 2,035 22,524 
Table 3-7. Acreage Бу crop/use type within the Lower West Coast Planning Area. 
Year Citrus Crops Nursery Pasture Sod Sugar Golf Landscape 
2000 168,790 74,414 2,263 20,256 811 126,550 13,297 14,714 
2004 157,863 70,558 2,380 18,472 360 131,345 16,902 21,792 
2010 142,645 74,185 5,517 23,190 527 136,140 19,655 29,620 
2014 125,798 73,267 5,116 22,556 681 126,289 20,924 32,843 


AG and L/R permittees often use several sources of water. Within the LWC Planning Area, permittees can 
withdraw from surface water and groundwater sources, including the Lower Hawthorn and Upper Floridan 
aquifers. Therefore, it was essential to create facility ratios that allowed for irrigation demands to be applied 
to the correct aquifers. Facility ratios were created by querying the SFWMD’s regulation database for 
historical water use data from 1999 to 2014, which were sporadic for AG and L/R permittees and had 
several outliers. A thorough quality assurance process examined each pumpage record and attempted to 
adjust any incorrect values. Following the quality assurance process, a facility ratio was created by taking 
the summation of each individual facility’s withdrawal and dividing by the total permit withdrawal. This 
allowed for each facility, regardless of the source of water, to have a ratio. 


Once the ratios were created, the SFWMD’s regulation database was queried for facility information, 
including x and y coordinates, pump capacity, and year in addition to well-specific information such as 
cased depth, total depth, and well diameter. Facilities were filtered to exclude wells using the Lower 
Hawthorn and Upper Floridan aquifers, which are not part of the LWCSIM. Also excluded were facilities 
pumping water from external sources (e.g., canals, lakes, reservoirs) or on-site surface water sources 
hydraulically connected to external sources. Finally, facilities were filtered to include only primary wells 
and pumps. This filter aligned with the assumption that permittees use only facilities that already had a 
pump installed. Facilities were placed in the model based on the x and y coordinates in the SFWMD’s water 
use database. If the x and y coordinates were missing or incorrect, the facility was placed in a reasonable 
location within the permit boundary. When historical pumping reports were not available, the facility ratio 
was assumed to be equally distributed among the primary facilities, as indicated in the SFWMD's water 
use database. 


AFSIRS Model for Estimating AG and L/R Water Use 


AFSIRS was developed for Florida's water management districts by the University of Florida Institute of 
Food and Agricultural Sciences as a method to determine agricultural water use allocations for consumptive 
use permitting programs (Smajstrla 1990). The model estimates irrigation requirements for Florida crops, 
soils, irrigation systems, and climate conditions. The irrigation requirement for crop production is the 
amount of water, excluding precipitation, that should be applied to meet a crop's supplemental demand 
requirements without a significant reduction in yield. 
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The AFSIRS model is based оп a water budget of the crop root zone and the concept that crop ЕТ can be 
estimated from potential ET and К.. The water budget includes inputs to the crop root zone from rain and 
irrigation, and losses from the root zone by drainage and ET. The water storage capacity in the crop root 
zone is the product of the soil’s water-holding capacity and the depth of the effective root zone for the crop 
being grown. Irrigation is scheduled based on an allowable level of soil water depletion from the crop root 
zone. This level of simulation model development produced a functional model capable of addressing the 
wide variety of crops, soils, and irrigation systems in Florida and including variations in daily rainfall and 
ET rates. 


Reference Crop and Crop Coefficients 


Reference crops were used to develop the water demands for various crops grown in Florida. Daily ET for 
each crop was calculated as the product of reference ET and the crop К, for that day. The K. values varied 
with crop growth stage. Crop К. values were obtained from the literature and revised by the Institute of 
Food and Agricultural Sciences for bahia grass, citrus, and sod. Three separate sets of К. values were 
examined and included in the AFSIRS updates for this modeling effort. 


АП individual and general irrigation permits issued within the LWCSIM domain were reviewed. 
Information obtained from the permits included crop type, irrigation method, irrigation source, and 
irrigation withdrawal facilities. The permitted average annual allocations also were collected and reviewed 
to ensure the AFSIRS-generated demand was less than the permitted annual allocation. Individual demands 
were calculated on a crop-by-crop basis and aggregated up to a monthly irrigation demand for each permit. 


For the transient model, estimated demands were applied to permitted groundwater withdrawal facilities. 
Surface water withdrawals were not included in the model unless the permittee used an isolated lake. In 
this case, the demand was applied to the SAS (layer 1 — Water Table aquifer). Irrigation demands were 
calculated from AFSIRS and varied based on climatic conditions. The SFWMD's water use database was 
reviewed to determine when a newly issued permit became operational within the simulation period. For 
these users, pumpage values before permit issuance were assigned a value of zero. 


3.3.4  Commercial/Industrial/Institutional 


CII water use is the smallest use category of the SAS/IAS within the LWCSIM domain, making up less 
than 10% of total water use. СП includes permitted commercial, industrial, and institutional users who are 
self-supplied and do not receive water from a PS utility. Examples of CII use types include concrete 
processing, mining, and food processing. Uses that return water back to the withdrawal source (e.g., rock 
washing, heating and cooling, equipment washing) are not included. Within the LWCSIM domain, there 
are 15 CII permits and 71 facilities (Figure 3-17). The largest CII users are the Ortona and Palmdale sand 
mines in Glades County. 


Within the LWCSIM domain, CII water use for the transient model was obtained on a well-by-well basis 
from reported pumping values submitted by permittees to the SFWMD or from estimates based on reported 
values. When only total permit pumping values were available, pumping was assumed to be equally 
distributed among all facilities associated with the permit. This allowed for creation of a master well file 
that provides individual pumping rates for each well. 
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Figure 3-17. Location of Commercial/Industrial/Institutional pumping facilities using the surficial and 
intermediate aquifer systems in the Lower West Coast Planning Area. 
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3.3.5 Domestic Self-Supply 


DSS wells in Lee, Collier, Glades, and Hendry counties and associated DSS well data were compiled from 
several sources. The SFWMD contacted county officials to request information on DSS wells, including 
location, casing depth, and total depth. Information was provided as a geographic information system (GIS) 
shapefile or an Excel spreadsheet that could be converted to a GIS shapefile, so all DSS wells could be 
located in the LWCSIM domain. Information from Lee and Collier counties was the most extensive and 
complete. DSS wells in Charlotte and DeSoto Counties (within the LWCSIM model domain) were obtained 
from SWFWMD, as discussed in the next section. 


АП DSS well data were combined into a comprehensive, consistent Excel table containing county, state 
plane coordinates, well casing depth, well total depth, and per capita use rate (PCUR). GIS was used to 
intersect these points with the LWCSIM grid to assign a row and column to each well. 


The PCUR was calculated from information in the 2012 LWC Water Supply Plan Update (SFWMD 2012). 
The planning document contains projections of DSS user populations and water demand for each county. 
Specifically, in Appendix A of the 2012 LWC Water Supply Plan Update, Tables A-1 through A-8 contain 
the data used to calculate a PCUR for all DSS wells in each county. According to the 2012 LWC Water 
Supply Plan Update (page 6 of Appendix A), the PCUR is the total annual water use divided by the number 
of permanent residents. 


For example, to calculate the net average PCUR for DSS wells in Lee County in 2005: 
1. DSS-using population = 68,566 permanent residents (SFWMD 2012, Table А-2) 


2. Net average DSS water demand = 8.37 million gallons per day (mgd) (SFWMD 2012, Table А-7) 
= 8,370,000 gallons per day (ера) 


3. PCUR = 8,370,000 gpd + 68,566 permanent residents = 122 gpd per person (gpd/p) 


The PCUR also was calculated for the net 1-in-10-year drought (129 gpd/p), the gross average (147 gpd/p) 
and the gross 1-in-10-year drought (156 gpd/p). The average of all the above is 139 gpd/p. 


Each county’s overall average PCUR based on the 2012 LWC Water Supply Plan Update (SFWMD 2012) 
was as follows: 


e Collier County: 219 gpd/p 

e Glades County: 139 gpd/p 

e Hendry County: 131 gpd/p 

e Lee County: 139 gpd/p 

Discussions with SFWMD planning staff indicated a better estimate of the PCUR might be obtained using 
the SFWMD’s Water Supply Utility Project (WaSUP) database of finished residential water demand as 
reported for each PS utility in the LWC Planning Area. Based on information in the WaSUP database, each 
county’s average PCUR was adjusted to the following: 


Collier County: 123.67 gpd/p 
Glades County: 114.58 gpd/p 
Hendry County: 67.43 gpd/p 
Lee County: 71.05 gpd/p 
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The derived PCURs were used as starting points for assigning extraction rates to model cells representing 
DSS wells. The rates per model layer per model cell were determined and assigned using a method 
developed by the SFWMD, which uses well casing depth and total depth to determine which model layer(s) 
are being tapped. The method accounted for model cells that contain more than one DSS well. It was 
decided that the PCUR represents water used per person and that each household well provides water for 
an average of three persons. Therefore, a factor of 3 was included when arriving at the final per-model-cell 
pumping rate. In other words, each DSS well represents three persons. Application to all four counties 
resulted in the following: 


e Collier County: 123.67 gpd/p x 3 = 371.01 gpd/well 

e Glades County: 114.58 gpd/p x 3 = 343.74 gpd/well 

e Hendry County: 67.43 gpd/p x 3 = 202.29 gpd/well 

e Lee County: 71.05 gpd/p x 3 = 213.15 gpd/well 

Based on information received from Lee County and discussions with SFWMD planning staff, two 
additional pieces of information were used to refine the DSS pumping rates. First, the DSS well average 
PCUR for each county was increased by a factor of one to account for private landscape irrigation, which 
resulted in the following: 


Collier County: 371.01 gpd/well + 123.67 gpd/well = 494.68 gpd/well 
Glades County: 343.74 gpd/well + 114.58 gpd/well = 458.32 gpd/well 
Hendry County: 202.29 gpd/well + 67.43 gpd/well = 269.72 gpd/well 
Lee County: 213.15 gpd/well + 71.05 gpd/well = 284.20 gpd/well 


Second, the pumping rates for all DSS wells were adjusted over time according to seasonal variations 
measured at DSS wells in a portion of Lee County. This information was obtained from a Lehigh Acres 
groundwater flow model developed by Copp (2015). Mr. Copp provided the SFWMD with calibrated DSS 
well pumping rates, partly based on private landscape irrigation use, and information on how these rates 
vary monthly over a 2-year period (January 2013 through December 2014). The monthly variation in DSS 
well pumping rates was used to derive a multiplier (ranging from 0.86 to 1.14) that was applied to the PCUR 
for each DSS well for each month of the LWCSIM calibration/verification period (January 1999 through 
December 2014). 


As a check on the above water extraction rate-per-well, the four-county average PCUR of 373 gpd/well 
divided by 3 persons/well = 124 gpd/p. This number is almost identical to the 2015 WaSUP finished 
residential water demand value for Collier County (123.67 gpd/p). 


The final PCUR per well per county was assigned to DSS wells in the LWCSIM cells, with adjustment for 
seasonal variations applied monthly from January 1999 through December 2014. In the LWCSIM domain, 
there are 59,340 DSS wells in Collier, Glades, Hendry, and Lee counties, represented by a total of 
22,139 model cells (Figure 3-18). Appendix C contains a detailed presentation of how the DSS coverage 
was developed. 
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Figure 3-18. Domestic Self-Supply well locations using the surficial and intermediate aquifer systems in 


the model domain. 
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3.3.6 Southwest Florida Water Management District Ритраде 


Different methodologies were used to prepare the Southwest Florida Water Management District 
(SWFWMD) pumpage data that are beyond the LWC Planning Area but still inside the LWCSIM domain. 
Therefore, these processes are discussed separately. 


Non-DSS Wells 


The SWFWMD maintains a pumpage database for the DWRM (Environmental Simulations, Inc. 2014). 
SFWMD staff received non-DSS (e.g., PS, AG, СІ) pumpage data from the SWFWMD (К. Peterson, 
SWEWMD, personal communication, January 2019). The data were plotted in ArcGIS, and wells within 
the LWCSIM domain were selected (Figures 3-15 to 3-17). Because the data were monthly, reported in 
cubic feet per day, and formatted for MODFLOW, no additional work was needed regarding the flows. 
However, the layering needed to be adjusted to be compatible with the LWCSIM. For DWRM, layers were 
assigned based on open interval and geographic location (Environmental Simulations, Inc. 2014). Nine 
geographic zones were used because several wells had poor well construction information: 


Zone 1: Upper Floridan aquifer (UFA) and Lower Floridan aquifer (LFA) 
Zone 2: UFA and LFA 

Zone 3: SAS, UFA, and LFA 

Zone 4: SAS, UFA, and LFA 

Zone 5: Permeable Zone 3 (PZ3) and UFA 

Zone 6: UFA and LFA 

Zone 7: UFA and LFA 

Zone 8: PZ2, PZ3, and UFA 

Zone 9: SAS, UFA, and LFA 


According to Figure 3-12 from Environmental Simulations, Inc. (2014), the SWFWMD portion of the 
model area falls into Zone 3. In addition, Figures 3-13 to 3-15 were used to evaluate the model layering. 
SFWMD staff attempted to follow the layering methodology as much as possible. However, issues that 
complicated this process included the following: 


e There is little geologic control in the area. 

e Some model layers thin or pinch out in the SWFWMD. 

e Some wells are completed in confining units, while the LWCSIM team agreed to assign 
withdrawals to aquifers only. 


The following strategy was used for layering: 


a) Wells open to a single aquifer layer were assigned to that layer. Most of these were in layers 1 or 
9, but a few were in layer 7. 

b) Wells open to layer 4 were reassigned to layer 1. Layers 2 and 3 are thin or absent in these areas. 
Thus, layer 1 is the most likely source. This is consistent with the DWRM. 

c) Wells open to multiple zones in the SAS were assigned to layer 1. This is consistent with the 
DWRM. 

d) Wells open to layers 5 through 7 were assigned to layer 7. This corresponds to layer PZ1 in the 
DWRM. 

e) Wells open to layers 7 through 9 were assigned to layer 9. This corresponds to layer PZ2 in the 
DWRM. 

f) Wells withdrawing below layer 9 were considered to be withdrawing from PZ3 and/or the UFA. 
These withdrawal zones are below the LWCSIM base. 
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DSS Wells 


The SWEWMD hired Janicki Environmental to prepare the DSS data. The following is an extraction from 
a memo from Janicki Environmental (2018) to SWFWMD staff, describing the process (using SAS 
program) for preparing DSS data for the DWRM: 


“1. Extract data from WMIS and subset based on established selection criteria 

Fix STR’s (section, township, and range) where necessary in historic database 

Calculate proportion of total county wells in each county STR 

Read in county usage and apportion STR usage based on number of wells in each STR 

Merge in STR centroid information 

Export as a txt file to read into GIS (Invariant DSS annual for GIS.txt) to final deliverables folder 
Apply seasonality coefficients to create monthly timeseries for each STR 

Export to a txt file (Invariant DSS monthly.txt) to final deliverables folder” 


SAAN рор 


A similar procedure was applied to the DSS wells to prepare them for the LWCSIM. An additional step 
was needed because the DSS data were given in million gallons. The data had to be converted to cubic feet 
per day. 


4 MODEL CONCEPTUALIZATION AND DESIGN 


4.1 Software Selection 


The MODFLOW code (Harbaugh 2005), which is an industry standard computer program for 
finite-difference, three-dimensional groundwater flow modeling, was used for this modeling effort. 
Specifically, MODFLOW 2005 containing the SFWMD wetland module, was used. MODFLOW was 
selected for the following reasons: 


It has been widely accepted in the groundwater modeling profession for more than 30 years. 

The code is well documented and within the public domain. 

The code is readily adaptable to a variety of groundwater flow systems. 

The code is modular and easily facilitates any modifications required to enable its application to 

the types of unique groundwater flow problems encountered in South Florida. 

5. It has been used extensively in existing groundwater flow models built for SFWMD water supply 
planning areas. 

6. The SFWMD version of MODFLOW 2005 has been customized to include the Wetlands package. 


4.2 Model Boundaries 


To [ue 


Boundary conditions can determine model flow patterns and the model's water budget. Therefore, defining 
the model boundaries is an important aspect of model design. In order to minimize potential errors during 
the simulation, physical or hydrologic features are selected as the model boundary conditions when 
available. When unavailable, specified head boundaries are selected a great distance away from areas of 
concern in the model domain. The primary imposed stresses in the LWCSIM domain are pumping wells 
(e.g., for AG or PS), and the model will be used to evaluate the impacts of current and future pumping from 
these locations. Figure 4-1 shows the locations of current pumping wells (excluding DSS wells) in the SAS 
and IAS. 


The western boundary coincides with the Gulf of Mexico and Peace River and includes Pine Island and 
Sanibel Island. The eastern boundary coincides with Lake Okeechobee and the SFWMD's primary water 
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management canals (from north to south, L-24, L-25, L-4, and L-28) out to Everglades National Park. The 
southern boundary coincides with the Lostmans River and Big Lostmans Bay tidal boundaries. The northern 
boundary is a specified head boundary located a reasonable distance from areas of concern in the model 
(about 20 miles away from the northernmost pumping well). 
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Figure 4-1. Location of water withdrawal facilities by aquifer in the Lower West Coast Planning Area. 
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4.3 | Spatial Discretization 


The LWCSIM domain covers an area from the southern edge of Highlands County at the north to Key 
McLaughlin in Monroe County at the south, and from Charlotte Harbor and the Gulf of Mexico on the west 
to approximately the western borders of Palm Beach and Broward Counties to the east. The model grid is 
aligned in a north-south orientation. The model grid is composed of 553 rows and 512 columns, with 
uniform grid spacing of 1,000 ft. Vertically, the model is composed of nine layers. There are 2,548,224 cells 
in the model. Of these, 1,670,609 cells are active and 877,615 cells are inactive or no-flow. Grid size 
selection was based on the planned use of the LWCSIM, data availability, апа computational 
considerations. Based on the North American Datum 1983 State Plane Florida East coordinates (FIPS 0901 
feet), the model origin at the southwest corner of the model is: X-direction: 218,436; Y-direction: 441,788. 


Model layer definitions and the corresponding hydrogeologic units are provided in Table 4-1. The aquifers 
and confining units are treated as separate layers in the model. Future revisions to the model, or the 
development of sub-regional models, may require additional sub-layers to obtain local spatial resolution or 
computational stability. 


Table 4-1. Model layers and corresponding hydrogeologic units. 


Model Layer Hydrogeologic Unit Abbreviation 
1 Water Table aquifer WTA 
2 Tamiami confining unit TCU 
3 Lower Tamiami aquifer LTA 
4 Upper Hawthorn confining unit НІ 
5 Sandstone aquifer clastic zone SSA-clastic 
6 Sandstone aquifer confining unit SC 
7 Sandstone aquifer carbonate zone SSA-carbonate 
8 Mid-Hawthorn confining unit H2 
9 Mid-Hawthorn aquifer MHA 


4.4 Temporal Discretization 


The major stresses to the LWCSIM are boundary conditions related to recharge from rainfall, ET, recharge 
in the northern portion of the model, wellfield withdrawals, and seasonal tidal patterns. The primary purpose 
of the LWCSIM is to address long-term (20 to 50 years) planning issues, where a longer stress period could 
be used. А secondary application may include development of a companion tool for water use permitting 
purposes, which would require shorter stress periods. Considering the scale and resolution of the LWCSIM, 
the model's intended use, and temporal data availability, especially withdrawal data, monthly stress periods 
were selected. 


4.5 Calibration Period 


Because the LWCSIM will be used to predict the availability of water resources for the next 20 years or 
more, the transient model had to be calibrated to extreme conditions (e.g., floods, droughts) to have greater 
confidence in the model's predictive capability. Based on historical rainfall data in the LWC Planning Area 
from 1965 to 2014 (described in Section 2.1), average annual rainfall is 53.6 inches. The probability 
exceedance curve developed for this period shows average annual rainfall during 1-in-10 year drought 
conditions is 44.8 inches, while the average annual rainfall during 1-in-100 year drought conditions is 
approximately 41 inches (Figure 4-2). Representing 1-in-10 year drought conditions, 42.7 inches of rainfall 
were recorded in the LWC Planning Area in 2007, while there was 41 inches of rainfall in 2000, which 
represented 1-in-100-year drought conditions. 


56 


80 


70 


Rainfall (їп.) 


0 10 20 30 40 50 60 70 80 90 100 
Probability of Exceedence 


Figure 4-2. Rainfall probability exceedance curve for the Lower West Coast Planning Area, 1965 to 
2014. 


In the LWC Planning Area, available monitor wells for calibration purposes increased from 30 wells in 
1998 to 106 wells in 1999. Considering extreme dry conditions (droughts in 2000, 2001, and 2007), extreme 
wet conditions (hurricanes in 2005), and data availability (Figure 4-3), the 16-year period from 
January 1999 through December 2014 was selected as the transient model period of record. As a result, this 
approach provided insight to the potential changes in hydrologic conditions to meet projected needs during 
extreme conditions. Due to time and computational constraints, a shorter transient calibration period 
(2008 to 2012) was chosen for use in the Parameter Estimation software (PEST; Doherty 2010) transient 
calibration. For the transient calibration, 2008 was a wet year and 2012 was a dry year. For manual 
calibration, 1999 through 2012 were used (Figure 4-4). Years 2013 and 2014 were used for model 
verification. 
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Figure 4-3. Annual distribution of groundwater monitor wells in the Lower West Coast Planning Area. 
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Figure 4-4. Annual average rainfall for the transient simulation period in the Lower West Coast 
Planning Area. 
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5 MODEL IMPLEMENTATION 
5.1 Boundary Conditions 


Figure 5-1 depicts the general head boundary (GHB) around the perimeter of the active model domain, and 
the no-flow boundary (inactive cells) that lie outside the active model domain. Whenever possible, 
hydrologic features were used as boundary conditions. Figure 5-2 shows the surface water stations 
associated with the GHB cells in the WTA and TCU. The station names in Figure 5-2 coincide with the 
station names in DBHYDRO. 
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Figure 5-1. Location of the general head boundary, layers 1 through 9. 
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Figure 5-2. Surface water stations used to develop the general head boundary in the Water Table 
aquifer. 
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5.1.1 General Head Boundaries 


According to McDonald and Harbaugh (1988), a GHB consists of a water source outside the modeled area 
that supplies or removes water to a model cell at a rate proportional to the head difference in the active cell 
and the head assigned to the external source by the user. The rate at which water is supplied to a cell is 
given by equation (6): 


Q =C x (Hs — Hc) (б) 


Where: 


Q is the flow rate (ft*/day) to or from the boundary cell 

C is the constant of proportionality for the boundary (conductance) 
H, is the average head at the source boundary 

Нс is the average head in the cell 


Along the eastern border for the WTA (layer 1), initial conductance values were obtained from the LECsR 
Model (Giddings et al. 2006). For the northern, southern, and western boundaries, the conductance values 
were calculated using equation (7): 


C=KxB (7) 


Where: 


С is the conductance value (ft?/day) 
K is the hydraulic conductivity of the cell (ft/day) 
B is the thickness of the layer (ft) 


Hydraulic conductivities were interpolated from field measured data. Initial conductance values for the 
remaining layers also were determined by equation (7). 


Conductance values were adjusted by PEST during the calibration process. The stages for the boundary 
heads vary horizontally and vertically. For the WTA and TCU, the boundary heads were based on monthly 
stages at the surface water stations, if available. If there were no adjacent surface water bodies, the boundary 
heads were based on topographic elevations and physiographic regions (White 1970). For the remaining 
layers, the boundary heads were based on the monthly estimated stage for the GHB cells. For the LWCSIM, 
the Gulf of Mexico was considered a physiographic region and divided into five parts to facilitate the PEST 
calibration. 


MODFLOW is not a density-dependent model; therefore, water levels along the coast were adjusted based 
on total dissolved solids concentrations. There is a significant difference in density between fresh water and 


salt water when the total dissolved solids concentrations increase (Table 5-1). 


Table 5-1. Density differences of water samples relative to total dissolved solids concentration. 


Temperature (°С) Concentration (mg/L) Density (kg/m?) 
23 50 997.606 
23 500 997.948 
23 5000 1001.345 


°C = degrees Celsius; kg/m? = kilograms per cubic meter; mg/L = milligrams рег liter. 
Note: Calculations based on http://www.csgnetwork.com/h2odenscalc.html. 
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According to Figure 5-1, there are active cells within the Gulf of Mexico and other saltwater bodies. These 
areas are tidally controlled. Data from the National Oceanic and Atmospheric Administration (2015) were 
used to obtain tidal stages. According to the National Oceanic and Atmospheric Administration (2015) 
there are two active tidal stations within the LWCSIM domain: Fort Myers and Naples (Figure 5-3). 
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Figure 5-3. National Oceanic and Atmospheric Administration tidal stations. 
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Figure 5-4 is a hydrograph of the monthly stages at the Fort Myers and Naples tidal stations. Analyses of 
the data indicate the two stations have similar stages. Therefore, Fort Myers was chosen to represent stages 
in the tidal areas. 
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Figure 5-4. Average monthly stages (in feet mean sea level) at the Fort Myers and Naples tidal stations. 


Because MODFLOW is not a density-dependent model, water levels along the coast needed to be converted 
to an equivalent freshwater head. The following procedure was used to estimate the equivalent freshwater 
heads. Beach and Hood (2003) used equations (8) through (11) from Lusczynski (1961) to calculate the 
equivalent freshwater head. In addition, water level data from observation wells with total dissolved solids 
concentrations greater than 5,000 milligrams per liter (mg/L) were corrected with the same procedure. 


htf = hpf + Zi (8) 
hpf = sg(Zwl — (21с — Dcsg)) (9) 
Zi = Zls — Desg (10) 

Therefore: 
htf = (sg(Zw — (215 — Desg))) + (215 — Desg) (11) 
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Where: 


htf is the total freshwater head in the well 

hp is the pressure head in the well 

hpf is the freshwater pressure head in the well 
Zi is the elevation at base of casing 

Zls is the land surface elevation 

Zmsl is sea level, elevation = 0 

Zwl is the measured water level elevation 

Dcsg is the casing depth, below land surface 
Sg is the specific gravity of water in the well 


5.2 


5.2.1 


Initial Conditions 


Initial Heads 


Starting heads were prepared for the steady-state model layers using the following methodology. 
January 1999 was used as the starting date. 


1. 


Where: 


Groundwater levels and surface water stages were collected and compiled from databases 
maintained by the SFWMD (DBHYDRO), United States Geological Survey (USGS), and 
SWFWMD. Groundwater levels were assigned to the appropriate aquifer layer based on well 
construction and hydrostratigraphic layer definitions. Surface water levels were applied to layer 1 
to augment the groundwater data. 


Topographic data and a water level contour map from the nearby DWRM (Environmental 
Simulations, Inc. 2014) were used to supplement layer 1. 


The Geostatistical Analyst (a commercially available geostatistical tool) was used to determine the 
best interpolation method and compare the results with the actual data. 


Once reasonable results were obtained, the starting heads for the confining units were calculated 
by averaging the two adjacent aquifers. For example, the starting heads for the TCU were calculated 
using equation (12): 


Нтси = (Awra + Hira) = 2 (12) 


Нтси is the water level in the TCU 
Нуртал is the water level in the WTA 
Н;та is the water level in the LTA 


Figures 5-5 through 5-13 provide the starting heads in January 1999 for each model layer. These data were 
used for the steady-state model. Once the steady-state model was calibrated, the resulting heads were used 
for the transient version of the LWCSIM. 
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Figure 5-5. Starting heads of the Water Table aquifer. 
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Figure 5-6. Starting heads of the Tamiami confining unit. 
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Figure 5-7. Starting heads of the Lower Tamiami aquifer. 
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Figure 5-8. Starting heads of the Upper Hawthorn confining unit. 


68 


| 
SARASOTA Н 
Ят DESOTO 


т 
n 
« 
Іш 
© 
£ 
= 
< 
в. 


Starting Heads 


Layer 5 (Sandstone Aquifer - 
Clastic Zone) 
Feet 


ШШ 0-5 
ГЕЙ 6 - 10 
B] 11-15 
Г] 16 - 20 
[С] 21-25 
Г] 26 - 30 
ШЫ 31-35 
ГГ 36 - 40 
[С] 41-46 
Г] 47-51 
ЕЗ) LWCSIM Model Extent 


LWCSIM Active Model 
E Boundary 


MIAMI-DADE | 


LLL РА 5 ЖАК 
\\ad.sfwmd.gov\dfsroot\GS\GSBiz\WS\LWCSIM\Report_Figures\mxd\StartingHeads_SAC_Lyr5.mxd 


Figure 5-9. Starting heads of the Sandstone aquifer, clastic zone. 
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Figure 5-10. Starting heads of the Sandstone aquifer confining unit. 
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Figure 5-11. Starting heads of the Sandstone aquifer, carbonate zone. 
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Figure 5-12. Starting heads of the Mid-Hawthorn confining unit. 
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Figure 5-13. Starting heads of the Mid-Hawthorn aquifer. 
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5.3 Rivers 


The primary and secondary rivers/canals data set of the SFWMD’s Arc Hydro Enhanced Database and the 
USGS’s National Hydrography Dataset were overlaid to evaluate the consistency of the data. The 
rivers/canals attributes were triangulated and accessed from: 


1. Caloosahatchee River Basin Integrated Surface Water — Groundwater Model (SFWMD 2004). The 
model uses SFWMD surveys and past study documentations of canal design and surveys. 


2. SFWMD river/canal surveys. 


3. BCB Integrated Hydrological Model (DHI Water and Environment 2002, Atkins 2011). The DHI 
Water and Environment (2002) study covers the Tidal Caloosahatchee River Basin, Estero Bay, 
the BCB, the freshwater Caloosahatchee (C-43) River, and the Picayune Strand Restoration Project. 


4. The DHI Water and Environment (2002) report claims data were acquired from United States Army 
Corps of Engineers (USACE) and SFWMD surveys as well as past studies. 


5. South Florida Regional Simulation Model canal input data documentation, with data from SFWMD 
and USACE surveys and studies regarding: 


a. Central and Southern Florida Flood Control Project (C&SF Project) primary canals for the 
Lower East Coast Planning Area and water conservation areas, including selected secondary 
canals for the Lower East Coast Planning Area 


b. C-44 Canal attributes from the WaSh (URS 2008) modeling study 
Everglades Agricultural Area MIKE-SHE modeling study canals 
d. South West Florida Feasibility Regional Simulation Model study canals 


6. The LWC SAS model (Marco Water Engineering, Inc. 2006), which considered the canal detail 
designs from the USACE, multi-version canal GIS coverages used in the Regional Simulation 
Model, and canal cross-section data used in the MIKE11 model for a previous SFWMD project in 
the Naples area. 


Rivers in the LWCSIM principally include the SFWMD's primary and secondary canals and other perennial 
rivers and streams (Figure 5-14). In addition to rivers and streams, surface water bodies that were not 
included in the Wetlands package were included in the River package. Thus, Lake Trafford was modeled 
as a river. The LWCSIM uses the MODFLOW River package to simulate river flows. The MODFLOW 
River package was used instead of the more advanced MODFLOW Stream Flow Routing package to keep 
model execution and preprocessing times manageable. 


In MODFLOW, a river is a head-dependent boundary condition. The MODFLOW River package requires 
riverbed conductance and river stage as inputs to estimate the river-aquifer leakance. Riverbed conductance 
(C) is calculated by equation (13): 

C —(KXLXW)—M (13) 
Where: 
K is the vertical hydraulic conductivity of the riverbed sediment 
W is the width of the riverbed 


L is the length of the river segment within a model cell 
M is the thickness of the riverbed sediment 
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Figure 5-14. Location of rivers simulated in the model. 
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Cross-sectional and hydraulic characteristics were assigned to each model cell where rivers were identified. 
This included assigning reaches to model grid cells and determining river bottom elevations and area within 
each cell. As stated earlier, the cross-sectional data were triangulated from multiple sources, whereas K and 
riverbed sediment thickness data, except for Lake Trafford, were obtained from the calibrated values of the 
Marco Water Engineering, Inc. (2006) model. The calibrated K and sediment thickness were used as initial 
values for the LWCSIM. Aerial photographs were used for triangulation and estimation of the interpolated 
canal top widths. The missing bottom elevations of primary or secondary canals were estimated from 
USACE (1957, 1961, 1962, 1969) design books or GIS coverages from other SFWMD projects. River 
lengths were obtained from the National Hydrography Dataset or the Arc Hydro Enhanced Database. For 
Lake Trafford, bathymetric data were obtained from the SFWMD, sediment thickness was derived from 
ART Engineering, LLC’s (2004) geophysical survey, and K was estimated based on the nearest riverbed 
values. 


River stage data were needed to estimate leakance in and out of the rivers to the aquifers. When gauge 
stations were available, observed stages were used. However, in many cases river tributaries are not 
monitored; thus, various methods were used for river stage estimation depending on the circumstance. 


e Method 1: River stage level was interpolated for model cells between two gauge stations. Linear 
interpolation was based on river flow direction(s) and calculating a slope between two gauge 
stations’ water levels and the distance of each grid node from the origin of one of the stations. 


e Method 2: For some river grid cells, especially tributaries of main rivers where the upstream 
monitoring station is absent, water levels from shallow groundwater wells near tributaries on the 
upstream side were used as surface water level approximations. Where tributary rivers are 
connected to gauged rivers, the interpolated data of the model cell at the junction of the center of 
the tributary and the main river were used as the downstream data source. 


е Method 3: Some tributaries in the Caloosahatchee River do not have surface water gauge stations 
or shallow groundwater wells. In these cases, riverbank elevations, riverbed elevations, and aerial 
imagery of the rivers were analyzed. The stage in these cells was estimated to be 4 ft below the 
riverbank elevation. 


e Method 4: For the Cape Coral area, where a complex network of canals exists and water level 
slopes seem unrealistic, the stage was estimated as the average from multiple stations. 


In addition to areas where gauge stations were not present, there also were instances when river stage data 
were missing. For this modeling effort, incomplete stage data at multiple monitoring stations were filled 
using the method described in Appendix D (Chebud and Melesse 2011). 


5.4 Drains 


Drains in the LWCSIM represent most of the SFWMD’s tertiary canals in the region (Figure 5-15). These 
are surface water routing features and are restricted to the top layer of the model. The SFWMD’s Wetlands 
package was used in the steady-state and transient models to represent wetlands in the top layer. 


In MODFLOW, a drain is a head-dependent boundary condition. In a drain cell, the flow of water out of 
the aquifer depends on the assigned head and the conductance term. The head is compared to the computed 
head in the aquifer for the cell containing the drain. If the aquifer head is higher than the drain head, then 
the drain removes water from the aquifer. If the drain head is higher than the aquifer head, then the drain is 
considered to be dry and no water is removed or added to the aquifer. The amount of water removed is 
based on the conductance calculated by equation (13). 
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To develop the LWCSIM Drain coverage, two GIS coverages were used: 1) the National Hydrography 
Dataset, medium resolution, and 2) the digital elevation model for South Florida, light detection and ranging 
(LIDAR) 10 ft, found in the SFWMD’s Geospatial Data Library. The procedure used to extract the drain 
coverage for use in MODFLOW is described in detail in Appendix E. 
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Figure 5-15. Location of drains used in the model. 
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5.5 Wetlands 


The 2000 land use/land cover map (Figure 2-12) of the LWC Planning Area shows a substantial portion of 
the top layer of the model area is covered in wetlands, especially in the southern half of the model area 
(Figure 5-16). The Wetlands package was developed in 1998 by the SFWMD and the Hydrological 
Research Center at Florida Atlantic University to fulfill the need for a sound, physically based 
representation of wetlands and surface water/groundwater interaction. The Wetlands package accounts for 
sheet flow, wetting and drying, ET, and the vertical and horizontal flux components of the wetland/aquifer 
interaction (Restrepo et al. 1998). 
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Figure 5-16. Wetlands coverage in the model. 
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Surface water flow is treated in two broad categories: sheet flow through dense vegetation, and channel 
flow through a slough network. The Wetlands package uses a semi-empirical Manning-type equation rather 
than Darcy’s Law to represent surface water movement through dense vegetation. An explanation of the 
mathematical formulations for the Wetlands package may be found in Restrepo et al. (1998). Wetland cells 
are always simulated in layer 1. This package gives the user the option of simulating the wetland as an 
individual layer or embedding it within the aquifer. 


Channel flow in wetlands can occur through sloughs, canals, airboat trails, and lakes, which act as 
preferential flow paths for surface water in wetlands. Levees also can influence flow patterns in wetlands 
by acting as barriers to surface water flow. Preferential flow paths created by channels or local flow barriers 
can be simulated by specifying different anisotropic ratios for each wetland cell. This feature of the 
Wetlands package required modification to the original MODFLOW Block Centered Flow package, which 
holds the anisotropy factor constant throughout the layer. The establishment of preferential flow paths 
requires the user to specify which cell faces in the modeled feature are perpendicular to flow (in the case of 
barriers) or toward the direction of flow (for channels). Cell faces are numbered by convention as shown in 
Figure 5-17. 


Figure 5-17. Convention for numbering cell faces in the MODFLOW Wetlands package. 


Multipliers were applied to the selected cell faces to increase (in the case of channels) or decrease (in the 
case of levees) the interblock conductivity. Figure 5-18 represents a levee (in orange) acting as a barrier to 
flow occurring from left to right across the graphic, by specifying cell faces 2 and 3. In this case, face 3 is 
designated as primary and face 2 is secondary. 


Use of the Wetlands package in MODFLOW 2005 required several input data files defining the wetland 
cell, Kh of the soil type in that cell, vertical anisotropy, capillary fringe, and Kadlec number. These are 
prepared separately from the MODFLOW model, and the inclusion of a call file and the above data files in 
the MODFLOW name file incorporate the Wetlands package into the model run. Appendix F presents a 
detailed description of the GIS procedure used to build the wetlands coverage. 
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Figure 5-18. Example of a levee acting as a barrier to flow in Ше MODFLOW Wetlands package. 


5.6 Landscape Irrigation Return Flow 


Landscape irrigation (LSI) occurs generally in the residential/commercial mixed land use areas served by 
PS utility systems. LSI is an important component of the water balance that accommodates a portion of the 
reclaimed water flow and is combined with rainfall, potable water, and other irrigation water to satisfy 
landscape irrigation needs. The combined irrigation water applied to the land satisfies plant water needs, 
evaporates or is transpired, runs off, or infiltrates into the ground and percolates to the water table to 
complete the water cycle. 


PS utility service areas are the areas defined to receive LSI. The service areas are aggregated at the county 
level because of the complexity of characterizing the distribution of reclaimed water in the utility service 
areas for all the utilities in the LWC Planning Area. There often are instances where utilities’ potable water, 
wastewater, and reclaimed water service areas overlap, and other instances where utilities serve potable, 
wastewater, and reclaimed water outside their service areas via agreements between utilities. Adequately 
capturing these complexities in the time frame available for this modeling effort was deemed infeasible. 
Therefore, distributing LSI water at the county level is considered a reasonable simplifying assumption to 
address this issue. 


Potable and reclaimed water are the two sources of LSI water associated with PS systems. The portion of 
potable water applied as LSI (i.e., outdoor water use) was estimated by subtracting indoor water use from 
total PS pumping. For return flow, the total PS pumping includes not only pumping from the SAS/IAS but 
also from the Floridan aquifer system. Indoor water use is assumed to be equal to wastewater flow. The 
annual potable water use for LSI during the calibration period in the LWCSIM is shown in Figure 5-19. 
Only a portion of reclaimed water generated is used for LSI purposes, which was considered in this analysis. 
Historical wastewater treatment flows and reclaimed water use data were obtained from the Florida 
Department of Environmental Protection for municipal wastewater systems. Historical reclaimed water LSI 
quantities were extracted based on the assumptions noted above. Reclaimed water used for LSI during the 
calibration period for Lee and Collier counties is shown in Figure 5-20. These calculations were compiled 
at the county level for input to the LWCSIM. There are no major potable or reclaimed water uses outside 
Lee and Collier counties; thus, LSI return flows were not calculated for Charlotte, Hendry, Glades, or 
Monroe counties. 
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Once the LSI rates were developed, they were added to the rainfall array and processed through the CN 
method to estimate recharge and runoff. Spatial distribution of LSI application rates from 1999 are shown 
in Figure 5-21. Development of the LSI data set for the LWCSIM is described in Appendix G. 
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Figure 5-19. Potable water used for landscape irrigation in the model. 
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Figure 5-20. Reclaimed water used for landscape irrigation in the model. 
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Figure 5-21. Landscape irrigation application rates used to calculate infiltration by curve number method 
for 1999. 
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5.7 ET-Recharge-Runoff Methodology 


The methodology used to develop ET and recharge to the SAS uses AFSIRS (Smajstrla 1990) with the 
NRCS CN method for partitioning rainfall and runoff (Restrepo and Giddings 1994). A FORTRAN 
program was written to call AFSIRS for different land use polygons to calculate daily ET and recharge 
requirements, which are translated into model cell values (ET-Recharge-Runoff program, Bandara 2018). 
The ET-Recharge-Runoff program uses time-dependent data such as rainfall, irrigation return flows, 
potential ET, land use, and crop type, as well as time-independent data such as drainage basins, soil types, 
irrigated fractions, and irrigation efficiencies. Input requirements and how each of the processes are 
connected in the ET-Recharge-Runoff program are shown in Figure 5-22. 
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Figure 5-22. Inputs and outputs for the ET-Recharge-Runoff program. 
5.7.1 Input Requirements 


GIS Coverage Input: The ET-Recharge-Runoff program reads spatially referenced data that overlay the 
model grid. GIS coverages were created for land use, soil type, topography, and climate. These are the main 
inputs to the program and were produced using an area-weighted value for each model cell. 


Model Grid Arrays: The main inputs to the program related to model-grid-referenced data include basin 
identification (ID) for water budget calibration purposes, basin ID for runoff CN, and cell ID for 
active/inactive designation. (Note: CNs are empirical parameters used to estimate runoff from rainfall based 
on land uses.) 


Reference Land Use Table: Look-up tables relate land use to the following: crop type, CN, growing 
season, percent vegetation, impervious area, switching (indicating if a polygon is irrigated or not), and water 
use type classification. Irrigated land use categories include citrus, sugar cane, vegetables, pasture, and sod. 
Non-irrigated land use categories include rangeland and upland forests. 
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Evapotranspiration Input: The District s ET data set (Brown 2013) was used to assign reference ЕТ 
values for each cell in the model domain. The potential ET was computed from the reference ET and Ke. 


Rainfall Input: The District s NEXRAD rainfall data set (Brown 2014) was used to assign a rainfall value 
to each cell in the model domain. 


Return Flow Input: LSI flows and other applications that introduce water into the unsaturated zone were 
added directly into the rainfall arrays. Septic tank return flows were added directly to recharge arrays. 


5.7.2  Evapotranspiration 


ET is the sum of vaporization by the combined process of evaporation and transpiration (Jensen et al. 1990, 
Doherty et al. 2010). ET is the largest loss in the water budget in Florida. The ET rate from a crop depends 
on the type of crop, stage of growth, soil moisture content, and amount of energy available. The reference 
crop ET was used to represent a baseline rate of ET for actively growing grass, which has an unlimited 
supply of soil water. Potential ET for other crops was calculated as the product of the reference ET and К.. 


The reference ET data came from the SFWMD’s updated ET data set (Brown 2013). The К. value depends 
on several factors, including growth stage, development of the crop canopy, and amount of rainfall. К. 
ranges between 0.4 and 1.2 for most agricultural crops and varies from month to month. Low K. values 
occur early in the growing seasons when vegetation canopies are incomplete, whereas high К, values occur 
during peak growth stages. The K. value increases following rain or irrigation due to increased evaporation 
from wet soil. AFSIRS provides tables of monthly crop coefficients for perennial and annual crops. Daily 
К, values were calculated within the AFSIRS program using linear interpolation and adjusted based on the 
daily rainfall. Because the AFSIRS program is designed to estimate AG demands, it does not provide К. 
values for non-irrigated crop types such as upland forests. In such cases, the К. values were determined by 
comparing the crop type in the land use GIS coverages with the AFSIRS crop description. During the 
calibration process, the assumed К. values were adjusted, if necessary, to better match observed data. Total 
crop ET was calculated as the sum of unsaturated zone ET and groundwater ET (equation 14). 


ET, = ETyz + ETaw (14) 


Unsaturated zone ET was calculated within the AFSIRS program. The AFSIRS program does not directly 
output unsaturated zone ET, so the source code was modified to properly account for and output unsaturated 
zone ET. Once unsaturated zone ET is calculated, maximum groundwater ET was calculated. The 
maximum groundwater ET was provided as an input to MODFLOW, which used the simulated depth to 
the water table and the crop root zone depth to calculate the actual groundwater ET. Actual groundwater 
ET varies linearly with water table depth, with maximum ET rate at ground level and zero at the root zone 
depth. Groundwater ET for a model cell, Үер gr is calculated as: 


Veet, ET = YR (ЕТС pi E ETyzi) * Аретті (15) 


Where: Apery, is the pervious area of the land use class i. Pervious areas were calculated using data from 
the National Land Cover Database. In well-irrigated areas, the crop ET requirement was assumed to be met 
by irrigation. Therefore, in well-irrigated areas, groundwater ET approaches zero. In open-water conditions 
where surface water and groundwater are connected, the ET requirement was assumed to be completely 
met by groundwater and, therefore, groundwater ET approaches potential ET. 
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5.7. Runoff Calculations 


The CN method developed by the NRCS, formerly known as the Soil Conservation Service (1989), was 
used to compute runoff. Net rainfall (infiltration and retention) was computed as rainfall, plus return flows, 
minus runoff; however, additional criteria were used to consider antecedent conditions and high rainfall 
conditions. 


CNs were selected according to the antecedent moisture contents as per the Soil Conservation Service 
(1972, Table 42) classification. According to this classification, a dry condition was when the total 5-day 
antecedent rainfall was less than 0.5 inches during the dormant season or less than 1.4 inches during the 
growing season. A wet condition was when the total 5-day antecedent rainfall was more than 1.1 inches 
during the dormant season or more than 2.1 inches during the growing season. The condition in between 
these two are defined as normal conditions. CNs were selected from antecedent moisture content I, I, or 
Ш categories for the dry, normal, and wet conditions, respectively. 


If dry conditions existed, an additional retention factor was added to the NRCS initial abstraction. The 
retention factor is a user input and was adjusted during calibration. This procedure considers antecedent 
conditions to limit runoff when there is potential for additional soil storage. 


The majority of southwestern Florida soils are sandy with very high infiltration capacity; therefore, runoff 
generation occurs mainly due to a saturation excess mechanism. Runoff was calculated for each model cell 
and accumulated by watershed. Runoff volume for each model cell (Ү,гџ тор) was calculated аз: 


Veelt,rof f == Уе Drof fi * Aroffi (16) 


Where: 


i is the land use class within the model cell 
Drog fi 15 the runoff depth 

A rof fi is the area of runoff 

nLU is the number of land uses within the cell 


For this modeling effort, nLU was limited to three predominant land use categories. The area of runoff was 
the area directly connected to drainage features that contributed to the basin outflows. The directly 
connected areas were calculated using Sutherland's (2000) equations. Directly connected areas convey 
runoff from an originating point to a destination such as a drain or a river. Runoff volumes are obtained by 
multiplying runoff depths by directly connected areas. 


In open-water conditions, runoff was assumed to be negligible during the dry season and high during the 
wet season. Thus, the model calculated runoff only during the wet season, using a higher CN. In closed 
basins, where runoff does not leave the basin and ridge physiographic province areas and where runoff 
potential is very small, runoff was assumed to be zero. 


Effective Rainfall 
Effective rainfall is the portion of rainfall available for plant ET. It is calculated as rainfall minus runoff 
and initial abstraction. For closed basins and in ridge areas, runoff was assumed to stay within the 


watershed. Thus, the total rainfall was effective after accounting for initial losses in closed basins and ridge 
areas. 
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Routed Runoff (Muskingum) 


The ET-Recharge-Runoff program was used to transform the direct runoff (from NRCS) to a routed flood 
hydrograph using the Muskingum hydrologic routing procedure (Cunge 1969). The direct runoff computed 
with the NRCS equation for all areas was accumulated by basin in preparation for the Muskingum 
hydrologic routing procedure. 


The Muskingum method is a hydrologic river routing technique based on the continuity equation and the 
linear reservoir relationship. The Muskingum method is adaptable to a simple mathematical solution with 
variables such as weighting factors for inflow and outflow rates, a storage time constant, and travel time 
constant. The inflow hydrograph was made up of the direct runoff from each basin. The result of the 
Muskingum transformation was a routed hydrograph curve with an equivalent volume to the inflow 
hydrograph curve. The Muskingum parameters were calibrated by correlating routed runoff with actual 
observed flow values. 


5.7.4 Recharge to the Saturated Zone 


AFSIRS simulates a daily water balance in the unsaturated zone. Originally developed to determine 
irrigation demands, AFSIRS was adapted to simulate the unsaturated zone water balance in non-irrigated 
areas as well. Data for AFSIRS come from crop, irrigation, and soil databases. The crop database includes 
crop type, monthly crop coefficients, total and irrigated crop root depths, and allowable soil water depletion 
by month. The irrigation database contains irrigation system type, application efficiency of the irrigation 
system, percent irrigated soil surface area, and the ET fraction assumed to be extracted from the irrigated 
portion of the crop root zone. The soils database includes available soil type, vertical hydraulic conductivity 
(Kv), and porosity. 


The AFSIRS model is based on a water budget of the crop root zone, and the concept of crop ET can be 
estimated from potential ET. The water budget includes inputs to the crop root zone from rain and irrigation 
and losses from initial abstraction, runoff, deep percolation, and root zone/unsaturated zone ET. The water 
balance equation within the crop root zone of a polygon can be written as: 


ASTO = Rain + IRR — RCH — Runoff — IA — UZET (17) 


Where (in inches): 


ASTO is the change in soil storage 

IRR is irrigation (In non-irrigated areas, IRR becomes zero.) 
RCH is the deep percolation or gross groundwater recharge 
JA is the initial abstraction 

UZET is the unsaturated zone ET 


In AFSIRS, the control volume is the available water content of the plant’s root zone. If incoming rainfall 
(or irrigation) exceeds this control volume, the excess volume (“drainage”) is discarded as the combined 
process of runoff and recharge to the saturated zone. For this modeling effort, runoff (from the CN 
calculations) was removed from the daily rainfall (1.е., effective rainfall) and used as an input to AFSIRS. 
This made the AFSIRS drainage output equal to the recharge in the saturated zone, which was passed to 
the MODFLOW model. 
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Recharge volume for a model cell (Иец,ксн) is calculated as: 


Veel, RCH =; уер DRNi * Aperv,i + Raines f * (Aimperv,i ж DCIA;) (18) 


Where: Aimperv,i is the impervious area of the land use polygon i. The second term in the equation 
represents the groundwater recharge that can occur within impervious areas (e.g., urban developments) 
where runoff is routed to on-site lakes or ponds. 


For non-irrigated areas, vegetation used available water in the unsaturated zone and, if the plant's water 
extraction depth was deep enough, supplemented the additional demands from the water table. AFSIRS 
was used to calculate the overall potential crop ET (reference ET x К.) and the amount of water extracted 
from the unsaturated zone. Unsaturated zone ET is based on the daily water balance of rainfall, potential 
crop ET, and available soil water storage. The potential saturated zone ET is equal to the potential crop ET 
minus the unsaturated zone ET. The potential saturated zone ET was passed to MODFLOW, which 
calculated actual saturated zone ET based on the distance from the ET surface to the water table and on the 
extinction depth for a specific crop type. 


In sufficiently irrigated areas, the calculations were the same except irrigation was applied when the soil 
moisture storage was depleted to a specified level (the maximum allowable depletion). Thus, the potential 
saturated zone ET approached zero for irrigated areas. Nonetheless, the potential saturated zone ET for 
irrigated areas was passed to MODFLOW. Both ET and recharge were assumed to be negligible in 
impervious areas due to runoff. 


Figure 5-23 describes the steps of the ET-Recharge-Runoff program. The program was successfully 
verified against field measured data within the LWCSIM domain. 
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Figure 5-23. Flow chart describing the ET-Recharge-Runoff program. 
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Figures 5-24 and 5-25 show the spatial distribution of simulated runoff for 2005 (wet year) апа 2007 (dry 
year), respectively. Using the ET-Recharge-Runoff pre-processor, simulated runoff generally was lowest 
in the southern portion of the LWCSIM domain, which is primarily wetland areas such as Everglades 
National Park and Big Cypress National Preserve. This is because the overland flow over wetland areas 
were routed using the Wetlands package and runoff was not removed from the ET-Recharge-Runoff 
pre-processor. Simulated runoff generally was highest in areas with the most impervious land coverage 
(e.g., urban areas in Lee and Collier counties) and in areas where the soil is poorly drained and the water 
table is shallow (e.g., along the Caloosahatchee River floodplain). 


Figures 5-26 and 5-27 show the spatial distribution of unsaturated zone ET developed using the 
ET-Recharge-Runoff program for 2005 (wet year) and 2007 (dry year), respectively. Unsaturated zone ET 
is higher when the thickness of the unsaturated zone is greater. Unsaturated zone ET is lowest when the 
thickness of the unsaturated zone is small to non-existent or in shallow water table conditions (e.g., river 
floodplains). Generally, the unsaturated zone ET is higher in irrigated areas because the crop ET 
requirement is mainly met by the water content in the root zone and by the irrigation water, which is applied 
to the unsaturated zone. 


Figures 5-28 and 5-29 show the spatial distribution of maximum ET developed using the 
ET-Recharge-Runoff program for 2005 (wet year) and 2007 (dry year), respectively. Maximum ET is an 
input to MODFLOW and is highest in open water bodies and shallow water table conditions, while it is 
lowest in deep water table conditions. Figures 5-30 and 5-31 show the spatial distribution of gross recharge 
developed using the ET-Recharge-Runoff program for 2005 (wet year) and 2007 (dry year), respectively. 
Gross recharge is an input to MODFLOW and is highest in highly drained soil conditions, deep water table 
conditions, closed basins, or open water bodies. It is lowest in poorly drained soil conditions or shallow 
water table conditions. 


Figures 5-32 and 5-33 show the spatial distribution of MODFLOW-simulated groundwater ET for 2005 
(wet year) and 2007 (dry year), respectively. MODFLOW-simulated ET depends on the maximum ET 
provided (Figures 5-28 and 5-29), simulated depth to the water table, and extinction depth of the plants. 
Figure 5-34 shows the spatial distribution of extinction depths used for 2005 (Shah et al. 2007). Generally, 
the simulated groundwater ET is highest under open water or shallow water table conditions such as rivers, 
lakes, wetlands, and river floodplains. The simulated groundwater ET is lowest in deep water table 
conditions or in urban areas due to lack of tree cover. Figures 5-35 and 5-36 show the spatial distribution 
of simulated net recharge for 2005 (wet year) and 2007 (dry year), respectively. Net recharge is the 
difference between gross recharge and the MODFLOW-simulated groundwater ET. Generally, net recharge 
is highest in highly drained soil conditions and deep water table conditions. Net recharge is lowest or 
negative (aquifer is discharging) in poorly drained soil conditions and shallow water table conditions such 
as in the Caloosahatchee River floodplain. 
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Figure 5-24. Spatial distribution of simulated runoff in 2005 (wet year). 
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Figure 5-25. Spatial distribution of simulated runoff in 2007 (dry year). 
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Figure 5-26. Spatial distribution of unsaturated zone evapotranspiration in 2005 (wet year). 
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Figure 5-27. Spatial distribution of unsaturated zone evapotranspiration in 2007 (dry year). 
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Figure 5-28. Maximum evapotranspiration in 2005 (wet year). 
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Figure 5-29. Maximum evapotranspiration in 2007 (dry year). 
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Figure 5-30. Gross recharge in 2005 (wet year). 
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Figure 5-31. Gross recharge in 2007 (dry year). 
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Figure 5-32. MODFLOW-simulated evapotranspiration in 2005 (wet year). 
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Figure 5-33. MODFLOW-simulated evapotranspiration in 2007 (dry year). 
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Figure 5-34. Extinction depths in 2005. 
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Figure 5-35. Simulated net recharge in 2005 (wet year). 
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Figure 5-36. Simulated net recharge in 2007 (dry year). 
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5.75 Land Use Impact on Evapotranspiration and Recharge 


Spatial distribution of ET and recharge is largely attributed to the land use categories across the LWCSIM 
domain. Generally, more ET occurs in wetland and open water areas, while less ET occurs in urban and 
built-up areas. Highest net recharge occurs in managed recharge areas, usually within urban and built-up 
land use types, and lowest net recharge occurs in wetland areas. In southwestern Florida wetlands, 
groundwater generally discharges to the surface. Variation in simulated ET, recharge, and associated 
parameters with land use types was studied as part of this modeling effort to ensure the simulated values 
demonstrated generally accepted trends. Figure 5-37 shows the variation of gross recharge, runoff, 
unsaturated zone ET, actual ET, total ET, and net recharge calculated from the ET-Recharge-Runoff 
program with different land use types. FLUCCS codes are used to describe the land use types (SFWMD 
2009). Table 5-2 lists the major FLUCCS codes used for this modeling effort. 


Table 5-2. Florida Land Use Cover Classification System (FLUCCS) descriptions. 


Land Use FLUCCS Code Description 
1000 Urban and Built-Up 
2000 Agriculture 
3000 Rangeland-Upland Nonforested 
4000 Upland Forest 
5000 Water 
6000 Wetlands 
7000 Barren Land 
8000 Transportation, Communications and Utilities 


Figure 5-37a shows the variation of gross recharge among various land use types. Gross recharge 
represents the applied recharge to MODFLOW before ET losses were accounted for. As expected, the 
highest average gross recharge rates were simulated in open water conditions (FLUCCS code 5000; 
50 inches/year) and wetland areas (FLUCCS code 6000; 43 inches/year). Due to the absence of an 
unsaturated zone, gross recharge is mostly equal to rainfall in open water and wetland conditions. The third 
highest average gross recharge rates (24 inches/year) were simulated in urban and built-up (FLUCCS 
code 1000). Most built-up areas (e.g., residential communities) have on-site lakes or ponds that capture and 
store local runoff, which eventually recharges the water table. Irrigated agricultural areas (FLUCCS 
code 2000) have the lowest average gross recharge rates (14 inches/year). The greater ET requirement in 
agricultural areas makes less water available to infiltrate to the water table, reducing net recharge. 


Figure 5-37b shows average ET-Recharge-Runoff model-generated runoff for each land use type. The 
highest runoff (18 inches/year) was generated in the transportation land use type (FLUCCS code 8000). 
The transportation land use type generally includes highly impervious areas such as highways, roads, and 
airports; thus, greater runoff is expected. The second highest average runoff (17 inches/year) occurred in 
the upland-non-forested category (FLUCCS code 3000). This FLUCCS code is used in conjunction with 
various land uses, including rural residential, reclaimed mining lands, abandoned mining lands, inactive 
land with street pattern, utility corridors, and parks and zoos, which contain a lot of impervious areas. The 
third highest average runoff (15 inches/year) occurred in the urban and built-up land use category (FLUCCS 
code 1000). This category also includes highly impervious areas such as cities, towns, villages, shopping 
centers, and industrial and commercial centers, which generate higher runoff. 
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Figure 5-37. Gross recharge, runoff, unsaturated zone evapotranspiration, groundwater 
evapotranspiration, total evapotranspiration, and net recharge variation with land use. 
(Refer to Table 5-2 for the FLUCCS codes.) 


Figure 5-37с shows average unsaturated zone ET calculated by the ET-Recharge-Runoff program for 
various land use types. For open water conditions (FLUCCS code 5000), unsaturated zone ET was zero 
because of the absence of an unsaturated zone. Unsaturated zone ET was smaller (2 inches/year) in the 
wetland land use type (FLUCCS code 6000), possible because the unsaturated zone is completely absent 
or present only during the dry season. Unsaturated zone ET was highest (31 inches/year) in irrigated 
agricultural areas (FLUCCS code 2000) because of the supplemental crop ET demand applied as irrigation 
to the unsaturated zone. Average actual groundwater ET simulated by MODFLOW for different land use 
types is shown in Figure 5-37d. As expected, the highest groundwater ET was simulated in open water 
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conditions (FLUCCS code 5000), followed by wetlands (FLUCCS code 6000). Irrigated agricultural areas 
have the smallest groundwater ET (2 inches/year) because crop demand was met by applied irrigation, 
rainfall, and unsaturated zone water content. Urban built-up (FLUCCS code 1000) and transportation 
(FLUCCS code 8000) land use types also had low groundwater ET (3 inches/year). This is mainly because 
both land use types have higher impervious percentages, which minimizes the area available for 
groundwater ET to occur. Additionally, most of the pervious portions are landscaped and receive irrigation 
return flows to meet the crop ET demands. Total simulated ET, which is the sum of unsaturated zone ET 
and groundwater ET, is shown in Figure 5-37e. The highest total ET (58 inches/year) occurred in open 
water conditions (FLUCCS code 5000), followed by the wetland land use type (FLUCCS code 6000; 
44 inches/year), primarily due to the higher groundwater ET contribution in these two categories. Generally, 
the average total ET simulated in all other land use types varied between 28 and 37 inches/year. Average 
simulated net recharge for various land use types is shown in Figure 5-37f. The highest average net recharge 
(20 inches/year) was simulated in the urban built-up (FLUCCS code 1000) areas. Net recharge was highest 
due to the combined effect of relatively higher recharge (managed recharge) and minimal groundwater ET 
in urban and built-up areas. Simulated net groundwater discharge was approximately 8 inches/year in open 
water conditions, where ET is highest and recharge is minimal because of the saturated conditions. The 
second lowest net recharge occurred in wetland areas (1 inch/year), also due to higher groundwater ET rates 
and lower recharge rates. 


5.7.6 Extinction Depths Variation with Land Use Types 


Extinction depth, or root zone depth, is a key parameter that controls ET. Extinction depths generally vary 
with crop type and thus land use type. Variation of average extinction depths with different land types are 
shown in Figure 5-38. The extinction depths in the model were highest (11 ft) for open water condition 
(FLUCCS code 5000). The highest extinction depths are assigned to open water to allow for maximum 
groundwater ET extraction in MODFLOW. The second highest average extinction depth (8 ft) was assigned 
to upland forest (FLUCCS code 4000) because roots generally grow deeper to extract groundwater. The 
lowest extinction depth (3 ft) was assigned to irrigated areas (FLUCCS code 2000), where crop water 
demand is mainly met with applied irrigation, in addition to the unsaturated zone water content and rainfall. 
Therefore, roots generally do not extend deeper to extract water from water table. The groundwater table is 
kept lower in most irrigated areas to avoid root damages due to saturation. The transportation land use type 
(FLUCCS code 8000) was assigned a smaller value (3 ft) because those areas are landscaped and well 
irrigated. Crop ET demands are met by the application of LSI water and roots generally do not extend 
deeper to extract water from the water table. 
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Figure 5-38. Average simulated extinction depths for different land use types. (Refer to Table 5-2 for 
the FLUCCS codes.) 


5.8 Pumping Wells 


Pumping wells were entered into the model using the MODFLOW WEL package. The WEL package was 
generated using well input tables for each water use type: PS, CII, DSS, and AG with L/R. PS, CII, and 
DSS well input tables were constructed from measured or estimated pumping data, as described in 
Section 3.3. 


Monthly reported PS and СП data were obtained from the SFWMD’s regulatory database. However, 
reported data for AG and L/R were incomplete. The AFSIRS program was used to generate the monthly 
AG and L/R data by permit. AFSIRS-generated demands were capped at 0.8 times the permitted allocation 
because the AFSIRS-generated demands generally are 1.2 to 1.5 (1 + 0.8 = 1.25) times higher than the 
permitted allocations. AFSIRS demands were generated for average climatic conditions for the simulation 
period, while the permitted allocations represent 1-in-10 drought condition demands. Therefore, AFSIRS 
demands should always be less than permitted allocations. During the calibration process, AFSIRS demands 
were further adjusted to better match observed water levels. Monthly permit-level demands were assigned 
to well demands based on historical ratios, if available. If historical records were not available, the 
permit-level pumpage rates were equally distributed among the number of active wells. 


A FORTRAN program was developed to read the individual well files and generate the MODFLOW WEL 
package. This program reads monthly pumpage data from each pumping well, total and casing depths of 
the wells, K data, and aquifer top and bottom elevations and assigns the withdrawals to the correct model 
layer. If a pumping well is open to multiple aquifers, the program assigns pumping rates proportional to the 
transmissivity of the aquifer. The use of well input tables facilitated easy editing of pumping rates, as 
needed, especially during manual calibration. 
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6 CALIBRATION 


Model calibration involved adjusting model input parameters within reasonable ranges until model results 
(simulated water levels and/or structure flows) closely matched observed field data, while staying within 
the error bands of hydraulic properties observed from APTs or other information. For the LWCSIM, 
calibration was achieved by adjusting the parameters both manually and automatically. Manual model 
calibration was the process of changing one or more parameters in the model for each calibration model 
run. In the automatic calibration process, the model parameters were adjusted and optimized using the 
computer software PEST (Doherty 2010, Doherty and Hunt 2010). 


Both the steady-state (1999 conditions) and the transient models were calibrated. Calibration and 
verification of the LWCSIM was intended to represent the hydrologic conditions from 1999 through 2012. 
This transient calibration period showed significant changes in population, land use, and water use as well 
as multiple drought and wet cycles. Prior to the calibration process, calibration goals were identified that 
describe reasonable tolerance limits for the goodness-of-fit of the simulation results to the measured field 
conditions. For this modeling effort, the comparisons were spatial and temporal. Multiple adjustments to 
aquifer hydraulic property types and values and to water recharge and discharge related inputs were made 
in a focused, trial-and-error process until the simulation results reasonably achieved calibration goals. The 
resulting calibrated model simulated historical and future aquifer conditions within the limits of calibration 
and model construction. Monthly groundwater level calibration targets were developed for the LWCSIM 
domain for the calibration and verification period from January 1999 through December 2014 using the 
District’s DBHYDRO database and external county sources. 


Calibration of the LWCSIM proceeded in an iterative fashion, using both automated (PEST) and manual 
adjustment of model parameters. The overall calibration process, from start to finish, consisted of the 
following steps: 


1. A steady-state model representing 1999 conditions was constructed, tested, and calibrated using 
PEST, varying only Kh in all nine model layers. Layer pinch-outs were not incorporated into the 
model for this first PEST calibration. 


2. The results of step 1 were used as initial conditions for a full transient PEST calibration of the 
model, varying K; river, drain, and GHB conductance; ET and recharge; and aquifer storativity. 
Layer pinch-outs were incorporated for the PEST transient calibration. PEST transient calibration 
was successful in calibrating the SAS layers (layers 1 through 3). 


3 Manual calibration picked up where PEST calibration left off and was used to complete calibration 
of the LWCSIM, including the IAS layers (layers 4 through 9). 


These were the major steps of the calibration process. of the following sections describe each calibration 
step, from the initial steady-state PEST calibration to the final manual transient calibration. 
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6.1 Methodology 


6.1.7 Calibration Targets 


The calibration targets were primarily groundwater monitor wells that had observed data. Calibrating a 
groundwater model only to groundwater levels does not guarantee a unique solution, especially when 
groundwater systems are highly connected to the surface water system, as they are in southwestern Florida 
where surface water bodies are strongly connected to the SAS through wetlands and rivers. This 
necessitated ensuring the surface water fluxes and water levels also were reasonable. An iterative 
groundwater recharge estimation process using the ET-Recharge-Runoff pre-processor and MODFLOW 
ensured the modeled surface water fluxes were reasonable compared to the observed structure flow data in 
terms of overall mass balance and the transient response. Measured structure flow from two major 
watersheds were selected as flow calibration targets (Figure 6-1). Measured water levels at 60 gauge 
stations in wetland areas were used as surface water level calibration targets (Figure 6-2). Additionally, 
four measured ET sites were selected as ET calibration targets (Figure 6-3). There were 441 groundwater 
monitor wells used as calibration targets, including 297 in the WTA, 72 in the LTA, 29 in the SSA-clastic, 
and 18 in the SSA-carbonate, and 25 in the MHA (Figure 6-4). The goal of model calibration was to match 
the model-calculated head levels, flows, and ET to the measured data at these target locations. Time-series 
data for groundwater heads and surface water stages were collected primarily from USGS, SFWMD, 
SWEWMD, Lee County, and Collier County monitor wells. Measured ET data were obtained mainly from 
the USGS. The monitor well data were thoroughly reviewed, and questionable well targets were removed. 
Observed data were summarized as averaged monthly heads. Estimated baseflows at the measured structure 
flow locations also were used as calibration targets. Boundary fluxes from the eastern and northern 
boundaries of the LWCSIM were compared with the District’s LECsR Model (Giddings et al. 2006) and 
the SWFWMD's DWRM (Environmental Simulations, Inc. 2014), respectively. 
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Figure 6-1. Flow calibration targets. 
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Figure 6-2. 
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6.1.2 Automated Calibration 


The LWCSIM used a combination of automated and manual calibration methods. PEST (Doherty and Hunt 
2010), a model-independent parameter optimization code, was used first for automated calibration. To 
minimize computational time during PEST simulations, BeoPEST (Hunt et al. 2010) was used. BeoPEST 
manages parallel runs more efficiently than the standard PEST calibration method. 


Uniqueness in model calibration is achieved through a simplification of parameters called regularization. 
For this model, regularization was achieved through application of Tikhonov constraints, which solves 
ill-posed inverse problems. Tikhonov regularization is the imposition of a “smoothing constraint” on 
parameters (Doherty 2010). This approach applies each parameter difference to PEST as extra 
“observations.” Regularization also applies expert knowledge to the parameter estimation processes. For 
the LWCSIM, this was accomplished by constraining the aquifer parameter values while considering APT 
values and prior information gained from previously developed models. 


Pilot points as well as watershed-based approaches were used in the PEST calibration. Pilot points were 
used for parameter point values such as K and storage coefficients. The point values of these parameters at 
calibration are used to interpolate throughout the model domain using kriging. The pilot point technique 
helped to characterize heterogeneities in the hydraulic properties. Pilot points were generated following the 
guidelines in Doherty (2010): 


1) Place pilot points at the center of head calibration targets 
2) Place pilot points at locations of APTs and core samples (for Kv analysis) in the confining units 
3) Fill gaps in the model domain with initial educated guess values 


4) Refine the overall pilot point distribution, especially in dense areas prescribed by steps 1 and 2, 
resulting in an excessive number of pilot points over a small region 


A logarithmic transformation was applied to parametric values such as K, as the values range over several 
orders of magnitude. Watershed-based parametric calibration was done for parameters such as river and 
drain conductance. 


For the LWCSIM calibration effort, the PEST objective function (phi) to be minimized was defined as the 
sum of the weighted root mean square error of the heads and recharge fluxes. Weighting factors for each 
individual observation value were determined based on data reliability, number of observations, accuracy, 
and calibration performance. 


6.2 Calibration Parameters 


6.2.1 Hydraulic Conductivity 


The LWCSIM was calibrated by adjusting Kh values within the major aquifers and Kv values within the 
confining units. Flow is dominated by Kh in a productive aquifer and Kv in the confining units between the 
major aquifers. Parameters selected for calibration in PEST included Kh in the permeable aquifer units 
(1.е., WTA, LTA, SSA-clastic, SSA-carbonate, and MHA), Kv in the confining units, and storativity in the 
confined aquifers. In the LWCSIM, most head targets and APTs are located in the WTA, LTA, SSA-clastic, 
and SSA-carbonate. Therefore, most pilot points were generated and distributed in these units. In the MHA, 
aquifer parameter and head observations become much sparser spatially across the domain and are primarily 
concentrated along the coast. Pilot points in the confining units were spatially distributed in a grid-based 
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fashion. During the PEST calibration process, model input parameters were adjusted within user-specified 
ranges. The ranges for pilot points in the permeable units were set at reasonable multiples of observed 
values based on data confidence. During iterative calibration, some pilot point intervals were slightly 
modified. The parameter ranges were set to represent a bracketed range from observed field test values. 


A preliminary test run of PEST was performed in Groundwater Vistas to evaluate calibration of the 1999 
steady-state flow model to observed water levels using only Kh as the variable parameter. The test run 
yielded results that suggested maintaining continuous Kh distributions across each model layer, as opposed 
to representing hydrostratigraphic layer absences by assigning properties from the layer below in areas 
where the layer is apparently missing. Attempts to capture local heterogeneities adversely affected 
steady-state model stability. Thus, the layer pinch-outs were not incorporated in the steady-state calibration. 
Transient model calibration began with a K distribution that did include layer pinch-outs. This is discussed 
in more detail in Section 6.4.1. 


6.2.2 Leakance Values 


Leakance values (Kv) of the four confining layers—layer 2: the TCU, layer 4: the Upper Hawthorn 
confining unit (H1), layer 6: the Sandstone aquifer confining unit (SC), and layer 8: the Mid-Hawthorn 
confining unit (H2)—were adjusted by including them as a variable in the automated calibration process. 
The calibration method developed for the LWCSIM modeling effort relied on a regional approach, where 
important model parameters such as Kh and Kv were adjusted in a systematic way to best match observed 
water levels. 


6.2.3 Storage Parameters 


Storage parameters for the aquifer systems were estimated as a continuous spatial field variable instead of 
a single parameter across a model layer. Storage parameters help describe the rate at which water is taken 
from the aquifer matrix and help develop the time-varying responses of the transient model. Changes in 
these parameters were accomplished by enabling estimation of the spatial variation, such as the approach 
used with the Kh values of the aquifer layers estimated in the automated calibration process. 


6.2.4 Crop Coefficients 


The ET-Recharge-Runoff program calculates maximum groundwater ET, recharge, and runoff as described 
in Section 5.7. The maximum groundwater ET was calculated as the supplemental crop water demand, after 
using the water from rain, irrigation, and soil moisture. The program used the potential ET of each crop 
type associated with the land use and К. to calculate the crop ET requirement. By varying the crop ЕТ 
requirement, maximum groundwater ET demands varied. 


6.2.5 Extinction Depths 


Extinction depth is the cutoff depth below which actual groundwater ET cannot occur. The actual 
groundwater ET rate is highest at the ET surface (generally ground level), gradually decreases as it goes 
deeper, and ceases at the extinction depth. Generally, this depth is specific to crop type. Initial extinction 
depth estimates were obtained from previous modeling studies in the area (Shah et al. 2007). Extinction 
depths were varied during the calibration process to adjust ET rates. 
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6.2.6 Curve Numbers 


The NRCS CN method was used to separate runoff from rainfall, as described in Section 5.7. CNs are 
associated with land use, and higher CN values indicate higher runoff potential, while lower CNs indicate 
lower runoff potential. CNs were adjusted during the structure flow calibration process. 


6.2.7 River and Drain Conductance 


River and drain conductance values were adjusted to achieve better simulated baseflows, as described in 
Section 5.7. River and drain conductance values were joined with specific watersheds and assigned a 
watershed ID value for the grids situated within the watershed. For calibration, a similar multiplying factor 
was applied to each river/drain conductance within a given watershed ID. This avoided point-scale 
calibration of the thousands of grid cells, and instead downscaled it to a scale of four watersheds. Initial 
values of river/drain conductance were based on the values estimated using equation (13) in Section 5.3. 


6.2.8 Pumping Rates 


Pumping rates for AG, L/R, and DSS wells were estimated due to incomplete or poor quality data. Pumping 
rates were adjusted during the calibration process after the main calibration parameter changes were 
exhausted, particularly in situations where monitor wells showed trends in water levels (e.g., a sharp or 
gradual increase or decrease in water levels). DSS pumping rates were adjusted in the MHA (layer 9) 
pumping wells around Cape Coral to simulate the observed sharp declines in water levels over time. AG 
pumping rates in the SSA (layers 5 and 7) were adjusted during calibration to obtain a better match of 
simulated and observed water levels. 


6.3 Steady-State Calibration 


Calibration of the LWCSIM 1999 steady-state model was a multi-step process involving Groundwater 
Vistas, GIS, PEST++, and some manual adjusting of parameters. Water levels in 1999 from 
177 groundwater observation wells were averaged across all nine model layers. The model was run using 
the K value based on an interpolation of APT values, prior to PEST calibration, to ensure convergence. The 
results of this initial steady-state model run indicated an absolute residual mean of 2.95 ft, with an observed 
head range of 60 ft. PEST was then used, varying only Kh in the nine model layers. The initial Kh, based 
on APTs, was assigned to pilot points arrayed in an approximate rectangular grid across the active model 
domain. During the initial PEST calibration run, the Kh in each model layer was changed to achieve a better 
fit between observed and modeled heads. The total number of water level observation points was increased 
to 266 (with some original points omitted), with a head range of 46.14 ft. The absolute residual mean was 
reduced to 1.73 ft, and the residual statistics and modeled versus observed heads are shown in Table 6-1. 
The observed versus modeled head is presented in Figure 6-5. The resulting head solution was used as the 
initial head array for the start of the transient model run. In addition, hydrostratigraphic layer absences 
(i.e., pinch-outs) were incorporated into the model. 
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Table 6-1. Residual statistics after parameter estimation (PEST) calibration of steady-state model. 


ВИСШЕ Overall Statistics 
Layer Absolute Residual Mean | No. of Observations 
1 1.50 186 Residual mean -0.46 
2 0.98 11 reas 2. deviation = 
solute residual mean 
3 232 2d Residual sum of squares 1.27e4003 
4 1.93 12 Root-mean-square error 2.19 
5 3.11 14 Minimum residual -6.84 
6 4.17 4 Maximum residual 4.86 
7 161 8 Range of observations 46.14 
8 2.69 2 Scaled residual standard deviation | 0.046 
: Scaled absolute mean 0.038 
9 4.93 2 Scaled root-mean-square 0.047 
Overall 1.73 266 Number of observations 266 
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Figure 6-5. Observed versus simulated steady-state heads in all model layers. 


6.4 Transient Calibration 


For the transient model calibration, a combination of automated (PEST) and manual procedures was used. 
The main calibration targets were the monthly averaged observed water levels in 502 monitor wells and 
60 wetland gauges and the monthly averaged measured structure flows from 2 watersheds. Additionally, 
ET rates measured from four stations were used as calibration targets. Several model parameters were varied 
throughout the transient calibration process, such as K; river, drain, and GHB conductance; pumping rate; 


CN; crop coefficient; ET and recharge rate; and aquifer storativity. 


The calibration period ran from 1999 through 2012, while 2013 through 2014 were used as the verification 
period. The extent of flooded and dry cells, depth-to-water-table information, and regional potentiometric 
surfaces were used qualitatively to guide the manual transient calibration process. Additionally, GHB fluxes 
from saltwater boundaries were adjusted to minimize saltwater influx as artificial freshwater influx to the 


model. 
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6.4.1 Transient Calibration Setup 


Prior to performing the transient model calibration, a K distribution based on the steady-state PEST 
calibrations, which included stratigraphic pinch-out areas in layers 1 through 7, was tested in the 1999 
steady-state model. For the first pre-PEST full transient model (168 stress periods), this K configuration 
provided the fastest and most stable convergence (approximately 45 seconds per stress period). 


In order to incorporate the layer pinch-outs in layers 1 through 7 K value surfaces, Groundwater Vistas and 
ArcMap (GIS) were used. The procedure was as follows: 


1. Beginning with the PEST-derived К distribution in each layer, the К value was accessed апа 
post-processed using Groundwater Vistas. From there, the K value of layer 8 was exported as a 
point shapefile. 


2. The shapefile was brought into ArcMap and selected by location against a prepared shapefile 
representing layer 7 pinch-outs. 


3. The selected cells were extracted and saved as a shapefile containing layer 8 K values in the shape 
of the layer 7 pinch-outs. 


4. In Groundwater Vistas, the new shapefile was imported to the K value of layer 7. K values in the 
layer 7 pinch-out cells were replaced with the K values of the corresponding layer 8 cells. 


5. In Groundwater Vistas, the updated K value of layer 7 was exported as a point shapefile to use 
against layer 6, repeating the procedure above beginning with step 2. 


6. This process was repeated for the remaining model layers. 


The “bullseye” patterns or anomalies in Appendix Н, Figures H-1 to H-5, is the result of using the process 
above to insert adjacent layer properties in areas of layer pinch-outs. 


6.4. Transient Calibration Procedure 
Groundwater Level Calibration 


Calibration consisted of manually adjusting K, leakance and specific storage, pumping rates, ET, and 
recharge to improve matches between simulated and observed heads, flows, and ET rates. Starting 
conditions and initial parameters were obtained from the steady-state calibrated model for 1999. A 
continuing emphasis during the calibration process was honoring the conceptualization of the system 
described earlier. Parameters were constrained within reasonable ranges based on prior knowledge of the 
hydrogeologic system. 


PEST calibration focused on reducing the residuals of groundwater levels and baseflows by adjusting K, 
leakance, and river and drain bed conductance. PEST successfully met the calibration goals in the SAS 
layers of the model (layers 1 and 3, the WTA and LTA, respectively). However, no significant 
improvements were made by PEST to layers 5, 7, and 9 (SSA and MHA, respectively). Layers 5, 7, and 9 
were manually calibrated by adjusting K, leakance, storativity, and pumping rate. The leakance between 
the productive aquifers (i.e., between the MHA and SSA, between the SSA and LTA) was the most sensitive 
parameter during manual calibration. As expected, pumping rate in the MHA was a very sensitive parameter 
because the MHA is not a very productive aquifer. More detailed descriptions about parameter sensitivities 
can be found in Section 7. A sharp decline (on the order of 20 to 30 ft) in MHA water level was observed 
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from 2000 to 2003 due to the introduction of a large number of small DSS wells, primarily in the Cape 
Coral area. No reported data for the DSS pumping wells were available; therefore, pumping rates were 
adjusted until an adequate match between the observed and simulated water levels was obtained. 


Simulated ET rates were adjusted in wetland areas in the southern portion of the model to match measured 
ET rates. The ET rates were increased by increasing the wetland crop coefficients from 1.1 to 1.2. The 
simulated water levels in the southwest portion of layer 1 consistently were overestimated prior to the ET 
adjustments. After the ET adjustments, the simulated water levels were in good agreement with metered 
data. 


Structure Flow Calibration 


Surface water is highly connected to the SAS through rivers, drains, and wetlands in southwestern Florida. 
Therefore, it is important to ensure the interaction between groundwater and surface water systems are 
simulated reasonably well and the overall mass balance is preserved between surface water and 
groundwater systems. Simulated transient surface water flows were compared with measured surface water 
flows at gauged structures. The simulated total structure flow out of a watershed is the sum of runoff and 
baseflow within the watershed. Runoff was simulated using the ET-Recharge-Runoff program, and the 
simulated baseflow was extracted from MODFLOW river and drain leakages, as described in Section 5.7. 
Тһе ET-Recharge-Runoff program is not coupled to MODFLOW; therefore, a process was developed to 
iteratively calibrate canal net flows between the two programs. Interaction between the surface water and 
groundwater systems was facilitated by a simple canal water budget program and the ET-Recharge-Runoff 
program. The simple canal water budget program neglects storage in the canals and is represented by 
equation 19: 


[Runoff + Baseflow]~ner rLow = Structure canal outflows + Pumpage – (19) 
Structure canal inflows — Return flows from lakes 


Where runoff is the routed runoff computed with the NRCS/Muskingum method, structure canal inflows 
and outflows are obtained from measured flow at gauges, and baseflow (Q) is calculated using equation 20 
from MODFLOW: 


Qhaseflow = ((K X L X W) + M) х (he, — h) (20) 


Where: 


Бнуег is the river stage 

h is the aquifer hydraulic head 

K is the hydraulic conductivity of riverbed material 
L is the length of the reach 

W is the width of the river 

M is the thickness of riverbed sediments 


The first step was to develop initial ET and recharge values for MODFLOW using the ET-Recharge-Runoff 
program. MODFLOW was executed using the initial ET and recharge values, and the baseflow components 
were extracted from the River and Drain package outputs. For the selected basins, an instantaneous 
hydrograph of runoff was derived by subtracting estimated baseflow from observed stream flow hydrograph 
data and compared with runoff calculated from the ET-Recharge-Runoff program. Baseflows were 
estimated from net flows using the Web-based Hydrologic Analytical Tool (WHAT) program 
(https://engineering. purdue.edu/mapserve/WHAT/). The WHAT program uses frequency analysis of flows 
to conduct baseflow separation from managed watersheds. The WHAT results were compared against other 
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methods used by the USGS such as the Baseflow Index, stream flow hydrograph separation (HYSEP), and 
Partitioning to estimate daily discharge. The result shows the WHAT program has a better performance 
over the Baseflow Index for watersheds (larger than 500 square miles). 


Runoff and baseflow components were adjusted to achieve a reasonable match with the observed stream 
flow hydrograph. Runoff calibration was achieved by adjusting the CN of the model cells within the basin. 
Where the drainage basin area was significantly larger, the outflow hydrograph peak can be less pronounced 
and takes more time to occur than the instantaneous hydrograph. In such situations, runoff was routed using 
the Muskingum method to correct the time lag and attenuate the peak of the output flood hydrograph. 
Baseflows were adjusted based on the estimated values by adjusting the river and drain bed conductance. 
The ET-Recharge-Runoff program was re-run using the adjusted CNs, and then ET and recharge values 
were passed to MODFLOW. MODFLOW was re-run with adjusted river/drain bed conductance values. 
This procedure was repeated until a reasonable match between the observed and simulated runoff and 
baseflow hydrographs was achieved. 


Due to the size of the model cells (1,000 ft x 1,000 ft), there can be several different land use types within 
a cell. The initial CN for a given cell was calculated based on the three predominant land uses in the cell. 
Because CNs depend on the antecedent moisture content of the soil, they were adjusted based on the daily 
rainfall events during the calibration period, as described in Section 5.7. The model domain underwent 
substantial land use changes during the calibration period. To capture the effects of land use change on ET 
and recharge, four land use maps were used: 2000, 2004, 2010, and 2014. CNs were adjusted to reflect land 
use changes over the calibration period. 


Within the LWCSIM domain, two major watersheds were selected for calibration based on the availability 
of continuously measured structure flow data. Gauged basins were delineated by combining multiple 
HUC12 watersheds to ensure the availability of inlet and outlet stream gauging stations. 


Surface Water Level Calibration 


The LWCSIM routes overland flow over wetland areas, where the Wetlands package is used. Thus, wetland 
area water levels simulated with the Wetlands package can be directly compared to observed water level 
data. Initial Kh values were assigned based on kriged data, using the APT and other field data, and adjusted 
during calibration. Initial Kadlec coefficients were used from the LECsR Model and adjusted during 
calibration. In southwestern Florida, most wetlands are hydraulically connected to the water table; 
therefore, wetlands respond highly to groundwater changes. 


6.4.3 Transient Calibration Criteria 


Statistical calibration goals for the LWCSIM were established based on comparing modeled and measured 
water levels in 502 groundwater observation wells, providing monthly groundwater levels for the period of 
record. Not all 502 observation wells have data for every month of the period of record; therefore, not all 
wells were used in every model stress period. Model performance statistics were calculated by finding the 
differences between the measured and modeled values (referred to as the residuals) for each stress period 
and then assessing the result against the metric via direct comparison or calculating statistics on the 
residuals. 


The calibration goals for groundwater observation wells in LWCSIM were: 
1. А mean residual error for the five aquifer layers of less than 1 ft. 


2. Ап absolute residual mean and root mean square error of less than 5 ft. 
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3. 50% of the mean absolute simulated head residuals for all wells in the five aquifer layers had to be 
within 2.5 ft of observed. 


4. 80% of the mean absolute simulated head residuals for all wells in the five aquifer layers had to be 
within 5 ft of observed. 


5. 60% of the observed versus simulated hydrographs in the five aquifer layers had to have a Pearson 
coefficient of correlation (Е?) value greater than 0.4. 


6. 60% of the observed versus simulated hydrographs in the five aquifer layers had to have a 
transformed Nash-Sutcliffe efficiency (TNSE) coefficient value of greater than 0.4. 


R? is defined as: 


К = gi sy S100? (21) 


Where: 


S; = the i? simulated head at the monitor well 
О; = the 1° observed head at the monitor well 


5; = the average simulated head at the monitor well 


O; = the average observed head at the monitor well 


The TNSE was determined for all aquifer layer transient groundwater calibration targets as an additional 
indicator of transient model fit. Through statistical comparison with the R? value, it was determined that a 
TNSE value of 0 or higher roughly corresponds to an R? value of 0.4 or higher, allowing for determination 
of a TNSE calibration criterion for the LWCSIM that is similar to the R2. The TNSE was determined using 
the following equation: 


n .— 7.2 
TNSE 21- OLD (22) 


Where: 


X; = the i? transformed observed head value at a particular observation well 
Х;-О-0 


X = the mean of the transformed simulated head values at the observation well (this should be zero) 
7; = the i transformed simulated head value at the same well location 

Zi= Si-S 

5 = the mean of the simulated head values at the observation well 

5; = the average simulated head at the monitor well 


Oi, 0: and 5; as previously defined 


To be more precise, the В? quantile for an R2of 0.4 is about 0.262 (i.e., approximately 26% of В? values are 
less than 0.4). The TNSE quantile for a TNSE of about 0.063 is 0.262 (1.е., approximately 26% of TNSE 
values are less than 0.063). The TNSE quantile for a TNSE of zero is about 0.252. Based on the above, 
three categories of model fit were created based on the TNSE: 


е TNSE < 0: Model quality is less than null model of observed mean 


e 0O<TNSE < 0.5: Low-quality fit 
• 0.5 < TNSE: Moderate- to high-quality fit 
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Based on the above, а TNSE cutoff of 0 was chosen. 


The calibration criteria were derived from the recently calibrated, peer-reviewed East Central Florida 
Transient Expanded (ECFTX) Model (Central Florida Water Initiative Hydrologic Assessment Team 
2020). A summary of the transient calibration criteria is provided at the end of this section. 


An alternative method for using calibration statistics to assess the fit of the model to observed water levels 
was to use the range in head values as the basis for establishing calibration goals. For example, requiring 
the absolute mean head error and the root mean square error to be less than 5% to 10% of the range in head. 
In other words, if the range in head in any given month is 100 ft, then the goal would be 5 to 10 ft. 


The same calibration criteria were used for surface water level calibration. The calibration criteria for 
structure flows were selected based on the criteria used in previous regional models (Central Florida Water 
Initiative Hydrologic Assessment Team 2020) and are as follows. 


1. Deviation of volume less than 15% 
2. Nash-Sutcliffe efficiency coefficient greater than 0.5 
3. В? greater than 0.5 


Deviation of volume (DV) is defined as: 


N .—0; 
Dy = 2551220 100 (23) 
іші Vl 
Nash-Sutcliffe efficiency (NSE) is defined as: 
t (01-5)? 
NSE =1- с (24) 
R? is defined as: 
Nee. S20; 0272 
2 b Qi (Si-Si)(0i-01)) (25) 


COEPGUSO*SI (О-О)? 
Where: 
S; = the i^ simulated flow value at the structure 
О; = the 1° observed flow value at the structure 


5; = the average simulated flow value at the structure 
О; = the average observed flow value at the structure 


Final Calibrated Model Properties 


Appendix H contains figures showing the end result of the manual transient calibration: property 
distributions in the LWCSIM layers. Field test measurements of Kh are overlaid to show the calibrated Kh 
distributions stayed within an order of magnitude of field measured values during the calibration process 
throughout the LWCSIM domain. Transmissivity, leakance, and storage coefficient distribution figures also 
are presented in Appendix H. 
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6.4.4 Model Fit/Statistics 
Groundwater Level Calibration 


The result of transient calibration is presented in the summary of calibration statistics (Table 6-2). The 
majority of the monitor wells are open to the WTA (297 wells) and LTA (72 wells). All water level 
calibration criteria were achieved for all aquifers of the model except mean absolute error (MAE) in the 
MHA. The calibration criteria for average mean error is +1 ft for all aquifers. In the WTA, the mean error 
was -0.27 ft, indicating a slight underprediction of water levels. In the LTA, the SSA-clastic and 
SSA-carbonate had mean errors of 0.58 ft, 0.52 ft, and 0.7 ft, respectively, indicating a slight overprediction 
of water levels. The average mean error in the MHA was -1.08 ft, indicating an underprediction of water 
levels. The calibration criteria for standard deviation, MAE, and root mean square error were 5 ft for all 
aquifers. These three criteria were met for all aquifers except the MAE criteria in the MHA. Minimum and 
maximum residuals generally are 10 ft or lower in all aquifers. More than 90% of monitor wells in all 
aquifers showed less than 5 ft of error. At least 59% of monitor wells showed an error less than 2.5 ft in all 
the aquifers. The calculated R? and TNSE values indicate a very good transient response of the simulated 
water levels. The calculated R? values show that at least 6596 of monitor wells had R? values greater than 
0.4. At least 72% of the monitor wells had TNSE values greater than zero in all aquifers. 


Table 6-2. Transient model calibration statistics of the target monitoring wells in the model. 


Aquifer 
Senones EM WTA LTA |SSA-clastic) _ SSA- MHA 
carbonate 

Mean error +1 ft -0.29 0.58 0.46 0.61 -1.08 
Error standard deviation <5 ft 2.31 3.11 3.77 2. 2.45 
Mean absolute error (MAE) «5 ft 1.99 2.82 3.74 3.26 2.95 
Root mean square error «5 ft 2.32 3.14 3.73 2.15 2.63 
Minimum residual -9.12 -4.12 -9.91 -3.28 -7.38 
Maximum residual 6.45 9.35 5.27 3.78 4.39 
Number of observation points 297 72 29 18 23 
Percentage with MAE < 2.5 ft >50% 71 60 29 78 78 
Percentage with MAE < 5.0 ft >80% 97 92 90 100 91 
Percentage with R? > 0.4 >60% 76 81 69 83 65 
Percentage with TNSE > 0 >60% 78 89 83 89 12 


LTA = Lower Tamiami aquifer, МНА = Mid-Hawthorn aquifer; В? = Pearson coefficient of correlation; SSA = Sandstone 
aquifer; TNSE = transformed Nash-Sutcliffe efficiency; WTA = Water Table aquifer. 

Green font indicates compliance with calibration criteria. 

Calibration period: 1999 to 2012. 


Appendix I contains all calibration target hydrographs and a summary table for all model layers. Spatial 
distribution of groundwater monitor wells and their average simulated errors (residuals) for the WTA, LTA, 
SSA-clastic, SSA-carbonate, and MHA are shown in Figures 6-6 to 6-10. Monitor wells shown in red 
represent where average simulated water levels are lower than average measured water levels, indicating 
an underprediction. Monitor wells shown in blue represent locations where average simulated water levels 
are higher than average measured water levels, indicating an overprediction. The size of the symbols is 
proportional to the average residual at that location. Figure 6-6 shows monitor wells in the WTA (layer 1) 
are located mainly along the Gulf coast where the WTA is heavily used (e.g., Naples, Bonita Springs, Fort 
Myers, Cape Coral). Overall, the model shows a very good distribution of residuals with no significant 
spatial bias in the WTA. However, there is a consistent, slight overestimation of simulated water levels in 
the Picayune Strand area of southern Collier County. 
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Spatial distribution of groundwater target wells and the average errors of simulated heads in 
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Figure 6-7. Spatial distribution of groundwater target wells and the average errors of simulated heads in 
the Lower Tamiami aquifer (layer 3). 
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Figure 6-8. Spatial distribution of groundwater target wells and the average errors of simulated heads in 
the Sandstone aquifer — clastic zone (layer 5). 
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Figure 6-9. Spatial distribution of groundwater target wells and the average errors of simulated heads in 
the Sandstone aquifer — carbonate zone (layer 7). 
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Figure 6-10. Spatial distribution of groundwater target wells and the average errors of simulated heads in 
the Mid-Hawthorn aquifer (layer 9). 
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As shown in Figures 6-7 to 6-9, the LTA and SSA-carbonate/clastic aquifers show very good spatial 
distribution of errors with no significant bias. Most MHA monitor wells are in the Cape Coral area, and no 
significant bias is observed (Figure 6-10). The two MHA monitor wells (DBKeys 1660 and 37311) in 
southern Collier County indicate water levels are greater than 40 ft above the mean sea level, which is an 
artesian condition. The model, however, underestimates these water levels by more than 5 ft, and it was not 
possible to move the simulated water levels higher with reasonable parameter adjustments to match the 
observed levels. Observed water levels in the UFA (adjusted for equivalent freshwater heads) are slightly 
higher (approximately 42.5 ft) than MHA levels in this area. According to Geddes et al. (2015), confinement 
between the MHA and the Lower Hawthorn aquifer in western Collier County is not very thick. Bennett 
(2004) noted the separation between the MHA and UFA is weak, and there is an apparent hydraulic 
connection between the two aquifers in Collier County. Analysis of the nearest cluster wells open to the 
MHA (ВІСҮ-М71) and UFA (BICY-MZ2) indicates a very good correlation between the water levels in 
the two aquifers (Figure 6-11). This further confirms that separation between the MHA and UFA is 
semi-confined, suggesting there can be upward flow from the UFA to the MHA in western Collier County. 
However, the LWCSIM assumes no hydraulic connection between the two aquifers and uses a no-flow 
boundary at the bottom of the model, which may not be a true representation of the area. Therefore, the 
LWCSIM may be unable to accurately simulate the groundwater system in southwestern Collier County, 
and users are cautioned against relying on the model predictions of the MHA in this area. Because of this 
conceptual limitation, BICY-MZ1 and BICY-MZ2 were removed when calculating the layer-specific 
calibration statistics for the MHA. 
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Figure 6-11. Comparison of observed water levels in cluster wells open to Mid-Hawthorn aquifer 
(BICY-MZ1) and Upper Floridan aquifer (BICY-MZ2). 


Regression plots of simulated versus observed water levels for the WTA, LTA, SSA-clastic, 
SSA-carbonate, and MHA are shown in Figures 6-12 to 6-16. Most simulated data show good agreement 
with the observed data. Histogram plots of simulated versus observed water levels for the WTA, LTA, 
SSA-clastic, SSA-carbonate, and MHA, including a quantitative distribution of errors in each aquifer, are 
shown in Figures 6-17 to 6-21. Appendix J contains figures of the spatial distribution of TNSE and R? 
values for each aquifer. 
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Figure 6-12. Mean simulated versus observed water levels in the Water Table aquifer. 
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Figure 6-13. Mean simulated versus observed water levels in the Lower Tamiami aquifer. 
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Figure 6-14. Mean simulated versus observed water levels in the Sandstone aquifer — clastic zone. 
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Figure 6-15. Mean simulated versus observed water levels in the Sandstone aquifer — carbonate zone. 
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Figure 6-16. Mean simulated versus observed water levels in the Mid-Hawthorn aquifer. 
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Figure 6-17. Histogram of simulated versus observed water level differences in the Water Table aquifer. 
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Figure 6-18. Histogram of simulated versus observed water level differences in the Lower Tamiami 
aquifer. 
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Figure 6-19. Histogram of simulated versus observed water level differences in the Sandstone aquifer — 
clastic zone. 
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Figure 6-20. Histogram of simulated versus observed water level differences in the Sandstone aquifer — 
carbonate zone. 
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Figure 6-21. Histogram of simulated versus observed water level differences in the Mid-Hawthorn 
aquifer. 
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Surface Water Flow Calibration 


Figures 6-22 and 6-23 compare observed and simulated structure flow data for two major watersheds in 
the LWCSIM domain: West Caloosahatchee and East Caloosahatchee, respectively. Each figure shows the 
estimated versus simulated runoff, estimated versus simulated baseflow, total structure flow, and 
cumulative flow. Simulated runoff was obtained by subtracting the estimated baseflow from the total 
observed stream flow. Total observed flow at a watershed outlet was calculated as: 


Total watershed outflow = Observed flow at the basin outlet — Observed flow at the basin inlet 
Total simulated flow was calculated as: 
Total simulated outflow = Runoff + Baseflow from rivers and drains 


As described in the Structure Flow Calibration section, deviation of volume, the Nash-Sutcliffe efficiency 
value, and R? were used to evaluate the goodness of fit. 
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Figure 6-22. Observed versus simulated structure flow in the West Caloosahatchee watershed. 
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Figure 6-23. Observed versus simulated structure flow in the East Caloosahatchee watershed. 


The total simulated structure flow for the West Caloosahatchee watershed (Figure 6-22c) indicates the 
model captures most of the peaks during the simulation period. However, there are a couple years when the 
model was not able to capture peak flows (2001 and 2013), especially when the peak flows were 
significantly higher than average. These peaks are associated with how the inflow (S-78) and outflow (S-79) 
structures were operated during high rainfall/hurricane events, which were not simulated by the LWCSIM. 
The mass balance error (10.45% of deviation of volume in Figure 6-22d) occurred primarily due to not 
capturing these peak flows. The model indicated a significant groundwater contribution to the West 
Caloosahatchee watershed river flow during the calibration period. On average, the groundwater baseflow 
contribution is approximately 97 cubic feet per second (cfs), which is about 12% of the total stream flow 
(789 cfs) within the West Caloosahatchee watershed (Figure 6-22b). 


The total simulated structure flow for the East Caloosahatchee watershed (Figure 6-23c) indicates the 
model captures peak flows about half of the time during the simulation period. This is because the complex 
operational schedule of the inflow and outflow structures (S-77 and S-78) was not incorporated into the 
LWCSIM. In some situations, total outflow becomes negative when the outflow structure, S-78, releases 
less flow than the inflow structure, S-77, within the model time step (1 month). The model simulated mass 
balance error was approximately 3.72%. The model indicates that, on average, groundwater baseflow 
contribution is approximately 27 cfs, which is about 7% of the total stream flow (374 cfs) within the East 
Caloosahatchee watershed (Figure 6-23b). 
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Figures 6-22 апа 6-23 indicate the simulated total structure flows are in good agreement with observed 
data, especially given the complexity of structure operations was not incorporated into the LWCSIM. 
Table 6-3 shows calibration statistics of the West Caloosahatchee and East Caloosahatchee watersheds. 
Indicated by the green font in the table, both watersheds met all three established calibration criteria 
(deviation of volume, Nash-Sutcliffe efficiency, and R?). 


Table 6-3. Transient structure flow calibration statistics for major watersheds in LWCSIM 


Statistic Criterion East Caloosahatchee West Caloosahatchee 
Deviation of volume «1596 3.72% 10.45% 
Nash-Sutcliffe efficiency >0.5 0.55 0.51 


Pearson coefficient of 


correlation (R2) >0.5 0.58 0.52 


Green font indicates compliance with calibration criteria. 
Calibration period: 1999 to 2012. 


Surface Water Level Calibration 


A summary of surface water calibration at 60 wetland gauge stations is provided in Table 6-4. As indicated 
by the green font in the table, the LWCSIM met the established calibration criteria. Average mean error of 
all gauge stations was -0.76 ft, indicating a slight underestimation. Of the gauge stations, 85% had less than 
2.5 ft mean absolute error, 97% had less than 5 ft mean absolute error, and 63% had 60% of R? and modified 
Nash-Sutcliffe efficiency, indicating a very good overall transient response of the model. 


Table 6-4. Transient wetland water level calibration statistics in the model. 


Statistic Criterion Feet or % 

Mean error <+1 ft -0.76 
Error standard deviation <5 ft 1.98 
Mean absolute error (MAE) <5 ft 1.62 
Root mean square error <5 ft 2:14 
Minimum residual -8.52 
Maximum residual 4.26 
Number of observation points 60 

Percentage with MAE « 2.5 ft >50% 85 
Percentage with MAE < 5.0 ft >80% 97 
Percentage with R? > 0.4 >60% 63 
Percentage with TNSE > 0 >60% 63 


В? = Pearson coefficient of correlation; TNSE = transformed Nash-Sutcliffe efficiency. 
Green font indicates compliance with calibration criteria. 
Calibration period: 1999 to 2012. 


Spatial distribution of wetland gauge stations and their average simulated errors are shown in Figure 6-24. 
Hydrographs for the simulated and measured water levels are shown in the Appendix I. Gauge stations 
shown in red represent where average simulated water levels are lower than average measured water levels, 
indicating an underprediction. Gauge stations shown in blue represent where average simulated water levels 
are higher than average measured water levels, indicating an overprediction. The size of the symbols is 
proportional to the average residual at that location. 
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Figure 6-24. Spatial distribution of wetland surface water gauges and the average errors of simulated 
water levels. 
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Figures 6-25 to 6-28 compare simulated and measured ET data at four stations. The four stations represent 
different land use types, including agricultural land use (Figures 6-25 and 6-28), marsh land (Figure 6-26), 
and Cypress swamp (Figure 6-27). Simulated ET was adjusted using crop coefficients to better match the 
measured data. Figures 6-25 to 6-28 show a good agreement between the final simulated and measured ET 
data. 
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Figure 6-25. Measured versus simulated total evapotranspiration at the Ave Maria agricultural site. 
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Figure 6-26. Measured versus simulated total evapotranspiration at the marsh site. 
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Figure 6-27. Measured versus simulated total evapotranspiration at the Cypress swamp site. 
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Figure 6-28. Measured versus simulated evapotranspiration at the SGGEWX agricultural site. 
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6.5 Model Verification 


Once a model is deemed calibrated, a verification period outside of the calibration period is selected to 
evaluate the model’s performance and determine if the simulation continues to meet the calibration metrics 
previously assigned. In this case, the 2-year monthly transient period from 2013 through 2014 was selected 
as the LWCSIM verification period. 


Table 6-5 provides a summary of the LWCSIM verification statistics. Within the model domain, there were 
242, 48, 23, 14, and 15 monitor wells within model layers 1, 3, 5, 7, and 9, respectively, to compare 
simulated and observed water levels during the verification period. Calibration criteria were achieved for 
layers 1, 3, 5, and 7, and were comparable to the calibrated model; however, some of the calibration criteria 
were not met in layer 9. Layer 9 had fewer targets (15 targets during the verification period as opposed to 
25 targets during the calibration period) for the verification period, and many of those targets had significant 
data gaps. In addition, increased reclaimed water use in the Cape Coral area may have offset some of the 
DSS withdrawals in the MHA, making some DSS pumping rate assumptions invalid for the verification 
period. Layer 9 did meet the R? criterion and showed a residual mean error of 0.45 ft. In general, it can be 
expected that some targets would not meet the criteria during the verification period because calibration 
parameters can no longer be adjusted as they were during calibration. That very few wells did not meet the 
calibration criteria during the verification period indicates the model is reasonably calibrated and can be 
used for its intended purposes. 


Table 6-5. Transient model verification statistics of the target monitor wells. 


Statistics Criteria чыш 

WTA LTA SSA-2 SSA-1 MHA 
Mean error +1 ft -0.05 0.68 -1.14 -0.17 0.62 
Error standard deviation <5 ft 2.33 3.1 3.8 1.88 6.57 
Mean absolute error (MAE) «5 ft 1.95 2.79 3.13 2.36 Sod 
Root mean square error «5 ft 2.33 3.14 3.89 1.82 6.36 
Minimum residual -9.12 -4.12 -9.91 -3.28 -7.38 
Maximum residual 6.45 9.35 5.27 3.78 4.39 
Number of observation points 297 72 29 18 23 
Percentage with MAE < 2.5 ft >50% Jl 56 70 86 29 
Percentage with MAE < 5.0 ft >80% 97 90 91 100 27 
Percentage with R? > 0.4 >60% 83 85 74 86 71 
Percentage with TNSE > 0 >60% 72 86 87 93 64 


R? = Pearson coefficient of correlation; TNSE = transformed Nash-Sutcliffe efficiency. 
Green font indicates compliance with calibration criteria. 
Verification period: 2013 through 2014. 
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6.6 Simulated Heads 


Figures 6-29 to 6-33 present the simulated potentiometric surfaces of the five LWCSIM aquifers. These 
are average head surfaces produced for each aquifer from the calibrated transient model after a full 
192-stress-period model run. Figure 6-29 shows water levels in the WTA, which represents regional 
groundwater flow patterns reasonably well. The regional groundwater high is near Polk City in Polk County 
(outside the LWCSIM domain), and groundwater flows radially outward in every direction. Southeasterly 
flow from the Polk City high enters the LWCSIM through the northern model boundary. The regional 
groundwater flow that enters the model domain continues to flow southwest and eventually discharge to 
the Gulf of Mexico. The highest water levels along the northern boundary are 60 to 70 ft National Geodetic 
Vertical Datum of 1929 (NGVD29) near the Highlands-Glades county line. Another regional high water 
level area in the WTA is the Immokalee area, south of the Caloosahatchee River. The highest water levels 
in this area are 30 to 40 ft NGVD29, and groundwater flows radially outwards in every direction from this 
high. Northerly flow ends up in the Caloosahatchee River, while flows to the west and southwest discharge 
to the Gulf of Mexico. Low water levels along the Caloosahatchee River indicate groundwater flow towards 
the river from the north and south. The simulated potentiometric surfaces in the LTA, SSA-clastic, and 
SSA-carbonate shown in are similar to regional patterns in the WTA (Figures 6-30 to 6-32). The lowest 
potentiometric heads are in the MHA in the Cape Coral area, where water levels are -70 to -80 ft NGVD29. 


Figures 6-34 to 6-37 present the potentiometric surface (head) differences between the LWCSIM aquifer 
layers. The potentiometric head difference maps can help understand the regional flow vertical gradients 
within the model domain. These were produced by subtracting the average potentiometric surfaces between 
aquifer layers. For example, Figure 6-34 presents the average head difference between the WTA and LTA. 
This was produced by subtracting the average LTA head from the WTA head. Most of the time, WTA 
heads are 5 to 15 ft higher than LTA heads where the LTA is productive (i.e., within Lee, Collier, and 
Hendry counties). Figure 6-35 shows the simulated potentiometric head difference between the LTA and 
SSA-clastic. The potentiometric heads in the LTA are 10 to 15 ft higher in Lee County, indicating a flow 
gradient from the LTA to the SSA-clastic. Figure 6-36 shows the simulated potentiometric head difference 
between the SSA-clastic and SSA-carbonate. There is a potentiometric head difference from -10 to +10 ft 
within Lee County, indicating both upward and downward flow gradients. In the Cape Coral area, the 
potentiometric head in the SSA-clastic is approximately 10 ft higher than the head in the SSA-carbonate. 
Figure 6-37 shows the simulated potentiometric head difference between the SSA-carbonate and the MHA. 
In Cape Coral and Fort Myers within Lee County, SSA-carbonate heads are 15 ft higher than in the MHA, 
indicating a flow gradient from the SSA-carbonate to the MHA. In southern Collier County and eastern 
Glades County, MHA heads are greater than SSA-carbonate heads by approximately 10 ft, indicating an 
upward flow gradient from the MHA to the SSA-carbonate. 
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Figure 6-29. Simulated water levels for the Water Table aquifer. 
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Figure 6-30. Simulated water levels for the Lower Tamiami aquifer. 
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Figure 6-31. Simulated water levels for the Sandstone aquifer — clastic zone. 
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Figure 6-32. Simulated water levels for the Sandstone aquifer — carbonate zone. 
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Figure 6-33. Simulated water levels for the Mid-Hawthorn aquifer. 
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Figure 6-34. Simulated potentiometric surface differences between the Water Table aquifer and the 
Lower Tamiami aquifer. 
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Figure 6-35. Simulated potentiometric surface differences between the Lower Tamiami aquifer and the 
Sandstone aquifer — clastic zone. 
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Figure 6-36. Simulated potentiometric surface differences between the Sandstone aquifer — clastic zone 
and carbonate zone. 
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Figure 6-37. Simulated potentiometric surface differences between the Sandstone aquifer — carbonate 
zone and the Mid-Hawthorn aquifer. 
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67 Water Budget 


Water budget analysis was done for surface water and groundwater budget components. The surface water 
budget was analyzed using the budget terms produced in the ET-Recharge-Runoff program, including 
rainfall, runoff, gross recharge, unsaturated zone ET, maximum groundwater ET provided to MODFLOW, 
actual groundwater ET simulated by MODFLOW, net recharge, and change in soil water content. Water 
use withdrawal from all aquifers was included for comparison purposes. The groundwater budget was 
analyzed using the budget terms produced during the MODFLOW model run, including GHBs, wells, 
drains, rivers, recharge/ET, change in storage, and top/bottom flux. 


6.7.1 Surface Water Budget 


Using the ET-Recharge-Runoff program, average simulated water budgets for each county are shown 
Figure 6-38. Charlotte, Glades, and Hendry counties predominantly consist of agricultural areas; Collier 
County primarily consists of wetlands and agricultural areas; and Lee County mainly consists of urban 
areas. Agricultural and landscape irrigation combined (IRR in Figure 6-38) was highest in Hendry County 
(5.5 inches/year) followed by Lee County (3.5 inches/year). Hendry County’s total irrigation was mostly 
due to agricultural irrigation, whereas it was mainly due to landscape and golf course irrigation in Lee 
County. The rest of the budget terms are compared with respect to rainfall values because of different 
rainfall values in each county. 


Runoff rate was highest in Charlotte County, approximately 26% of rainfall (14.5 + 55 x 100 = 26%), while 
the lowest percentage of runoff occurred in Collier County, 15.7% of rainfall (8.9 + 56.4 x 100 = 15.7%). 
Although, the predominant land use is agriculture in Charlotte County, the portion within the LWCSIM is 
mostly urban. Therefore, runoff generated within Charlotte County is higher than the other counties within 
the model domain. Most of Collier County is covered with wetland areas. Runoff in wetland areas is not 
accounted for here because runoff is not separate from rainfall in the ET-Recharge-Runoff program; instead, 
runoff is routed through the Wetlands package. 


Unsaturated zone ET typically is highest in irrigated areas because it is assumed that the supplemental crop 
ET requirement is met by applied irrigation. Percentage of unsaturated zone ET with respect to rainfall was 
highest in Glades (25 + 48.8 x 100 = 51.2%), Hendry (23.3 + 52.7 x 100 = 44.2%), and Charlotte 
(21.1 +55 x 100 = 38.2%) counties. Groundwater ET rates were lowest in irrigated areas such as Glades 
(8.2 + 48.8 x 100 = 16.8%), Hendry (13 + 52.7 x 100 = 24.7%), and Charlotte (14.7 + 55 x 100 = 26.7%) 
counties. Groundwater ET was highest in wetland-dominated areas such as Collier County, 72% of rainfall 
(40.6 + 56.4 x 100 = 72%). Thus, total ET (unsaturated zone ET + groundwater ET) was highest in 
wetland-dominated areas such as Collier County, 83% of rainfall (46.8 + 56.4 x 100 = 83%), because of 
the higher contribution from groundwater ET. 


Net recharge rate was highest in Lee County, approximately 18.2% of rainfall (11.3 + 56.7 x 100 = 19.9%). 
Lee County also has the highest estimated water use withdrawal rates from the SAS and IAS, approximately 
3.4 inches/year. However, because of the greater net recharge rate (11.3 inches/year), water use withdrawals 
are unlikely to cause a significant threat to the SAS or IAS. Net recharge was lowest in Collier County, 
approximately 3.7% of rainfall (2.1 + 56.4 x 100 = 3.7%) because throughout much of the county, the water 
table exists either very close to or above ground surface, especially in wetland areas, causing groundwater 
to discharge to the surface. The estimated water use withdrawal rate from the SAS and IAS in Collier 
County was 1.8 inches/year, which was greater than the county’s net recharge. However, it could be 
misleading to use net recharge to represent the entire county because about half of Collier County is wetland 
areas. Therefore, comparison of water use withdrawals and net recharge may not be accurate. 
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Figure 6-38. Average water budgets, by county, as simulated in the ET-Recharge-Runoff program. 


6.7.2 Groundwater Water Budget 


A simulated average water budget, by model layer, from 1999 to 2014 is provided in Table 6-6. Average 
net recharge to the WTA was approximately 5.81 inches/year across the model domain (gross recharge — 
ET = 30.294 — 24.48 = 5.81). Average GHB lateral fluxes were 1.47 inches/year into the model domain. 
GHB lateral fluxes represent regional groundwater inflows to the model domain from the north/northeast 
minus groundwater discharge to the Gulf of Mexico. Average well withdrawal rate from the model domain 
was 1.842 inches/year. Most withdrawals came from the WTA and LTA. Total drain and river leakages 
from the model domain were 4.97 and 0.75 inches/year, respectively, mostly from the WTA. Drain and 
river leakages represent groundwater discharge as baseflows. Total average baseflow was approximately 
5.7 inches/year (4.97 + 0.720 + 0.025 = 5.715), which is comparable to the net recharge (5.81 inches/year). 
Generally, the SAS is highly connected to surface water bodies in southwestern Florida, so higher baseflows 
are expected. Total storage change was approximately 0.28 inches/year, indicating an increase in 
groundwater storage. This indicates groundwater withdrawals have not caused adverse impact to the 
SAS-IAS groundwater system overall. 
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Table 6-6. Simulated water budget summary. 


HB я Віуег n hange in Vertical Flow 
Layer s Wells | Drains m ro eund ET nodus Top RM Error 
1 0.573 | -0.378 | -4.974 -0.720 30.294 -24.480 0.274 0.000 -0.590 0.000 
2 0.000 | 0.000 | 0.000 -0.025 0.001 0.000 0.590 -0.566 0.000 
3 0.472 | -0.793 | 0.000 0.000 0.001 0.566 -0.245 0.000 
4 0.016 | 0.000 | 0.000 0.000 0.000 0.245 -0.261 0.000 
5 0.029 | -0.343 | 0.000 0.000 0.000 0.261 0.053 0.000 
6 0.019 | 0.000 | 0.000 0.000 0.000 -0.053 0.034 0.000 
Jd 0.017 | -0.237 | 0.000 0.000 0.001 -0.034 0.253 0.000 
8 0.048 | 0.000 | 0.000 0.000 0.000 -0.253 0.205 0.000 
9 0.297 | -0.093 | 0.000 0.000 0.000 -0.205 0.000 0.000 
Total | 1.470 | -1.842 | -4.974 -0.745 30.295 -24.480 0.276 1.117 -1.117 0.000 


ET = evapotranspiration; GHB = general head boundary. 
Note: All values are presented in inches per year. 


During the simulation period, the WTA provided an average recharge of approximately 0.59 inches/year to 
the underlying TCU, which eventually recharges the LTA. The WTA received 0.57 inches/year of net 
average groundwater recharge from the regional system from north/northeast of the model domain through 
GHB fluxes. Higher water levels along the northern/northeastern boundary and lower levels in the 
southwestern portion of the model indicate a northeast to southwest regional flow gradient (Figure 6-29). 
The WTA lost approximately 0.38 inches/year due to water use withdrawals, 4.98 inches/year through 
drains as baseflows, and 0.72 inches/year through rivers as baseflows. During the simulation period, there 
was a net increase in groundwater storage in the WTA of approximately 0.27 inches/year, indicating no 
negative impacts to the WTA from water use withdrawals. 


The LTA received 0.57 inches/year of recharge from the overlying WTA. Generally, the average 
potentiometric levels of the WTA are 5 to 15 ft higher than the LTA potentiometric levels where the LTA 
is productive (Figure 6-34). Therefore, a downward flow gradient is expected from the WTA to the LTA. 
The LTA received approximately 0.47 inches/year of net groundwater recharge from the regional system 
through the  northern/northeastern model boundary. Higher potentiometric heads along the 
northern/northeastern boundary and lower potentiometric heads in the southwestern portion of the model 
indicate a northeast to southwest regional flow gradient (Figure 6-30). The water use withdrawal rate from 
the LTA was approximately 0.378 inches/year, making the LTA the most used aquifer in southwestern 
Florida. Despite higher withdrawals, there was a slight increase in groundwater storage (0.001 inches/year) 
within the LTA, indicating no significant harm to the LTA from water use withdrawals. 


The LTA provided approximately 0.26 inches/year average recharge to the underlying SSA-clastic zone 
through the confining zone. Figure 6-35 shows the L TA potentiometric heads are 10 to 15 ft higher than 
the SSA-clastic zone potentiometric heads; therefore, a downward flow gradient from the LTA to the 
SSA-clastic zone is expected. The SSA-clastic zone received 0.03 inches/year of groundwater flow from 
the regional system, and 0.343 inches/year were withdrawn for water use purposes. No change in 
groundwater storage within the SSA-clastic zone occurred. The SSA-carbonate zone received 
approximately 0.032 inches/year of average recharge from the overlaying SSA-clastic zone and 
0.2 inches/year from the underlying MHA. The simulated potentiometric surfaces indicate there are areas 
in MHA with higher potentiometric levels than the overlying SSA-carbonate zone, especially in the 
southwestern portion of the model domain (Figure 6-37). Thus, an upward gradient from the MHA to the 
SSA-carbonate zone is expected. Approximately 0.24 inches/year of average water use withdrawals 
occurred from the SSA-carbonate zone; however, groundwater storage slightly increased 
(0.001 inches/year). 
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The MHA provided approximately 0.2 inches/year of average recharge to the overlying SSA. The MHA 
received 0.297 inches/year of net average groundwater recharge from the regional system through GHB 
fluxes from north and northeast of the model domain. Regional flow from northeast to southwest is clearly 
shown in the simulated potentiometric heads (Figure 6-33). Approximately 0.9 inches/year of average 
water use withdrawals occurred from the MHA; however, there were no changes in groundwater storage. 


Figures 6-39 to 6-43 present the average GHB fluxes into and out of the LWCSIM for the five aquifer 
layers. These represent one of the extracted MODFLOW budget terms and indicate the magnitude of 
groundwater flux entering and leaving the model under a regional groundwater gradient. Shades of green 
represent water flowing into the model domain and shades of yellow to red represent water flowing out of 
the model domain. In layers 1, 3, 5 and 7, the general flow across the model domain is from the northeast 
to the southwest, or from the Lake Okeechobee down and out to the Gulf of Mexico. This represents the 
known pattern of regional groundwater flow across the model domain. The GHB flux in layer 9 is mostly 
inflow from the Gulf of Mexico. There is a known area of artesian conditions in the southwestern portion 
of the MHA, as indicated by measured water levels in the two monitor wells in southwestern Collier County. 
Artesian conditions were simulated by defining higher stages in the GHBs along the Gulf coast, particularly 
in this area of the model. GHB boundary fluxes were checked to ensure there were no major fluxes on the 
western boundary indicating saltwater coming into the model domain, with the exception of layer 9. 
Western boundary fluxes show very little flow to the east because the corrected GHB heads for fresh water 
are higher than adjacent freshwater heads. 
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Figure 6-39. Average net general head boundary flow in the Water Table aquifer. 
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Figure 6-40. Average net general head boundary flow in the Lower Tamiami aquifer. 
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Figure 6-41. Average net general head boundary flow in the Sandstone aquifer — clastic zone. 
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Figure 6-42. Average net general head boundary flow in the Sandstone aquifer — carbonate zone. 
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Figure 6-43. Average net general head boundary flow in the Mid-Hawthorn aquifer. 
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7 SENSITIVITY ANALYSIS 


A sensitivity analysis is a process of quantifying the degree to which model outputs are influenced by model 
inputs. Model inputs include parameters that are conceptually constant or fixed, but whose values are not 
known with great precision. A sensitivity analysis may reveal the degree of calibration of a calibrated 
model. Under this concept, a sensitivity analysis is interpreted as a kind of inverse model; rather than 
estimating model performance for given inputs, error bars about the inputs are estimated. All models have 
some error. Useful models have errors that are not significant relative to the answers to the questions that 
are asked of the model. 


Whereas calibration is a minimization problem in which parameters are adjusted until model error is 
minimized, sensitivity analysis is descriptive; it explores the parameter space in the neighborhood of the 
calibration. Whereas the result of calibration is a single point, a list of specific parameters and their values, 
the result of sensitivity analysis is a collection of such lists. 


A model highly sensitive to a given parameter exhibits significantly increased error with deviation from the 
calibrated parameter value. Conversely, a model that is insensitive to a given parameter will not exhibit 
significant increases in error with such deviation. 


7.1 Parameters and Methodology 


Model sensitivity to several parameters was assessed by constructing several alternative parameter sets and 
computing relative model response. The result was several dozen pairs of deviations of parameter and 
response. Scenarios differ from the calibrated model by a single parameter, so the effect of that parameter’s 
variations on the model could be individually assessed. Parameter ranges were varied, within acceptable 
ranges, based on a predetermined data range for each parameter calculated from known values. Table 7-1 
shows the model input parameters used in the sensitivity analysis and the different multipliers applied to 
the calibration parameters. 


Table 7-1. Parameters varied by model layer. 


Horizontal hydraulic conductivity 


Each Layer* 
M Specific yield 


Each Pair of Adjacent Layers Vertical hydraulic conductivity 


Horizontal hydraulic conductivity of muck/peat 


Specific yield of muck/peat 


Specific yield of surface water body 


Kadlec conductance coefficient 


Layer 1 Drain conductance 


River conductance 


Recharge 


Evapotranspiration 


Evapotranspiration extinction depth 


Each of Layers 2 to 9 Specific storage 


* The distinction between “Each Layer" and “АП Layers” is that for the former, each layer is considered to have its own 
parameter, whereas for the latter, a single parameter is shared by all the layers. 
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Additional model runs replacing time-varying heads in selected GHB cells with a constant-head boundary 
in addition to a no-flow boundary (GHB conductance equal to 0) in model layers 1 through 9 also were 
simulated. For each parameter, several model runs were completed using the various multipliers identified 
in Table 7-2. The simulation period for the sensitivity analysis was from January 1, 1999 to December 31, 
2012 (168 months). To complete the sensitivity analysis, 178 model runs were conducted; of these, 
170 converged (ran to completion). 


Table 7-2. Parameters and multipliers of sensitivity analysis for the model. 


Description Multipliers 
Horizontal hydraulic conductivity 0.1 0.2 5 10 
(Layer 1) specific yield 0.25 0.5 0.75 1.25 
(Layers 2, 3, 4, б, 8, and 9) specific storage 0.01 0.1 10 100 
(Layers 5 апа 7) specific storage 0.0015 0.015 0.15 0.2 
(Layers 2, 3, 4, 6, 8, and 9) specific yield 0.05 0.1 0.2 0.25 
(Layers 5 and 7) specific yield 0.0015 0.015 0.15 0.2 
(Layers | to 8) vertical hydraulic conductivity * 0.1 0.2 5 10 
Drain conductance 0.01 0.1 10 100 
Evapotranspiration 0.8 0.9 1.1 1.2 
Evapotranspiration extinction depth 0.8 0.9 1.1 1.2 
General head boundary conductance 0.01 0.1 10 100 
General head boundary conductance > N=0 W= 0 E=0 i 
Recharge 0.8 0.9 1.1 1.2 
River conductance 0.01 0.1 10 100 
Kadlec conductance coefficient 4 0.1 0.2 5 10 
Horizontal hydraulic conductivity of muck/peat 4 0.1 0.2 5 10 
Specific yield of muck/peat ¢ 0.75 0.875 1.125 1.25 
Specific yield of surface water body 4 0.894 0.947 1.042 Ё 
а The model parameter varied for vertical hydraulic conductivity was actually leakance, but because the variation is by a factor, 
this is equivalent to varying vertical hydraulic conductivity. 


Three scenarios where the north, east, or west general head boundary conductance was set to zero. 
* Only three variations. 
4 Wetlands only. 


7.2 Sensitivity Results 


Sensitivity results are grouped by aquifer and subdivided into four sensitivity categories. Categories are 
defined by the magnitude of the difference between mean absolute error for the scenario and for the 
calibrated model (expression EXP1). 


| Sensitivity MAE – Calibrated Model МАЕ | (EXPI) 
Specifically, the difference categories are: 
е Less than 0.1 ft: not sensitive 
e From 0.1 to less than 0.5 ft: slightly sensitive 


e From 0.5 to less than 1.0 ft: moderately sensitive 
e 1.0 ft or greater: sensitive 
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A complete set of sensitivity plots, grouped by parameter type, is presented in Appendix К. There is one 
plot per parameter. The parameter value is plotted on the horizontal axis, and the responses (deviations 
from the base MAE for each aquifer) are plotted on the vertical axis. 


There is only one scenario for which there is an overall improvement in model performance: where Kh for 
layer 6 is reduced by a factor of five. However, this comes at the expense of an increase of relative MAE 
in layer | that approaches 2 ft. Table 7-3 provides additional details. In Table 7-3, the values are rounded 
to the nearest 0.1 ft. Layer 1 exhibits significant diminishment in performance for some stations. 


Table 7-3. Range of relative mean absolute error, by aquifer, for the scenario Kh (layer 6) = 0.2x. 


haven Range of Relative Mean Absolute Error 
Minimum Maximum 
1 -1.4 1.8 
3 -0.3 0.5 
5 -0.7 0.3 
7 -0.8 0.2 
9 -0.0 0.1 


Kh = horizontal hydraulic conductivity. 


7.3 Monitor Well Water Level Sensitivity 
Monitor well water levels were most sensitive to the parameters presented in Table 7-4. The LWCSIM was 
not sensitive in any aquifer to the scenarios not listed in the table. The model appeared most sensitive to 


Kh, Kv, and riverbed conductance. 


Table 7-4. Parameters with significant sensitivity. 


Parameter Layer | Multiplier цш Еш Rank 
WTA | LTA |SSA-clastic | SSA-carbonate | MHA Sum 
Kv 2 0.1 1 3 4 4 1 13 
River conductance N/A 0.01 2 2 2 2 4 12 
Kv 1 0.1 2 2 3 3 1 11 
Kh 3 10 2 3 2 2 2 11 
Kv 4 0.2 1 1 4 4 1 11 
Kv 4 0.1 1 1 4 4 1 11 
River conductance N/A 0.1 2 1 2 2 4 11 
Kh 1 0.2 2 2 2 2 2 10 
Kh 1 0.1 2 2 2 2 2 10 
Kv 2 0.2 1 2 4 2 1 10 
Kv 4 10 1 1 3 4 1 10 
Kv 4 5 1 1 3 4 1 10 
Kv 6 0.1 1 1 2 2 4 10 
Evapotranspiration N/A 0.8 2 2 2 2 2 10 
Kh 1 5 2 2 1 2 2 9 
Specific yield 1 0.25 2 2 2 1 2 9 
Kh 5 10 2 1 3 2 1 9 
Kh 5 0.1 1 1 4 2 1 9 
Kv 6 0.2 1 1 2 1 4 9 
Kh 9 10 1 1 1 2 4 9 
Kh 9 5 1 1 1 2 4 9 
Kh 9 0.1 1 1 1 2 4 9 
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Parameter Layer | Multiplier quite ЕЛІК К 

WTA | ГТА |SSA-clastic | SSA-carbonate | MHA Sum 
Drain conductance N/A 0.01 1 2 2 3 1 9 
River conductance N/A 10 1 2 2 2 2 9 
Recharge N/A 1.2 2 2 2 2 1 9 
Kv 1 0.2 1 2 2 2 1 8 
Kh 3 0.1 1 2 2 2 1 8 
Kv 5 0.2 1 1 1 1 4 8 
Kv 5 0.1 1 1 1 1 4 8 
Kh 6 10 2 1 2 2 1 8 
Kh 6 5 2 1 2 2 1 8 
Kv 6 10 1 1 1 2 3 8 
Kh 7 10 1 1 2 3 1 8 
Ку 7 10 1 1 1 1 4 8 
Ку 7 5 1 1 1 1 4 8 
Ку 7 0.2 1 1 1 1 4 8 
Kv f 0.1 1 1 1 1 4 8 
Kv 8 10 1 1 1 1 4 8 
Kv 8 5 1 1 1 1 4 8 
Kv 8 0.2 1 1 1 1 4 8 
Kv 8 0.1 1 1 1 1 4 8 
Kh 9 0.2 1 1 1 1 4 8 
Specific storage 9 100 1 1 1 1 4 8 
GHB conductance N/A 0.1 1 1 1 1 4 8 
GHB conductance N/A 0.01 1 1 1 1 4 8 
West GHB conductance N/A N/A 1 1 1 1 4 8 
River conductance N/A 100 1 2 1 2 2 8 
Specific yield 1 0.5 2 1 2 1 1 7 
Kh 2 10 1 1 1 2 2 7 
Kh 3 5 2 2 1 1 1 7 
Kh 3 0.2 1 2 1 2 1 7 
Kh 5 5 1 1 2 2 1 7 
Kh 5 0.2 1 1 3 1 1 7 
Kv 5 10 1 1 1 1 3 7 
Kv 5 5 1 1 1 1 3 7 
Kv 6 Э 1 1 1 2 2 7 
Kh 7 Э 1 1 2 2 1 7 
Kh 7 0.2 1 1 3 1 1 7 
Kh 7 0.1 1 1 3 1 1 7 
Drain conductance N/A 0.1 2 1 1 2 1 7 
GHB conductance N/A 100 1 1 1 1 3 7 
Recharge N/A 0.8 2 1 1 1 2 7 
Ку 2 10 1 2 1 1 1 6 
Ку 2 5 1 2 1 1 1 6 
Kh 4 10 1 1 1 1 2; 6 
Specific storage 8 100 1 1 1 2 1 6 
Specific storage 9 0.1 1 1 1 1 2 6 
GHB conductance N/A 10 1 1 1 1 2 6 
North GHB conductance | М/А N/A 1 1 1 1 2 6 
Evapotranspiration N/A 1.2 1 1 1 1 2 6 
Recharge N/A 0.9 1 1 1 1 2 6 


GHB = general head boundary; Kh = horizontal hydraulic conductivity; Kv = vertical hydraulic conductivity; LTA = Lower 
Tamiami aquifer; MHA = Mid-Hawthorn aquifer; N/A = not applicable; SSA = Sandstone aquifer; WTA = Water Table aquifer. 
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Rank was assigned by sorting the four sensitivity categories, with “not sensitive” ranked as 1. Based оп 
rank sum, among the parameters included in the sensitivity analysis, the LWCSIM was most sensitive to 
Kv for layer 2. A multiplier of 0.1 versus the base scenario effects an average of more than 1 ft in deviation 
from the base scenario (rank 4), in terms of MAE, for the SSA-clastic and SSA-carbonate. The LTA was 
moderately sensitive (rank 3). At least one model layer was sensitive (rank 4) or moderately sensitive (rank 


3) to Kv for layers 1, 4, 5, 6, 7, and 8. 


7.4 Results by Layer Parameter 


This section describes model sensitivity by aquifer, which is summarized in Table 7-5. Tables 7-6 to 7-10 


provide further details for each aquifer. 


Table 7-5. Summary of sensitive parameters by aquifer. 


Aquifer 


Sensitive Parameters 


Water Table aquifer (layer 1) 


Slightly sensitive to hydraulic conductivity, drain-river leakance, 
and recharge 


Lower Tamiami aquifer (layer 3) 


Moderately sensitive to vertical hydraulic conductivity between 
layers 3 and 4 and to its own horizontal hydraulic conductivity 


Sandstone aquifer — clastic zone (layer 5) 


Sensitive to vertical hydraulic conductivity between layers 2 
and 3, between layers 4 and 5, and between layers 5 and 6. 
Moderately sensitive (rank 3) to vertical hydraulic conductivity 
between layers 1 and 2 and between layers 4 and 5 

Moderately sensitive to horizontal hydraulic conductivity of 
layers 5 and 7 


Sandstone aquifer — carbonate zone (layer 7) 


Sensitive (rank 4) to vertical hydraulic conductivity between 
layers 2 and 3 and between layers 4 and 5 

Moderately sensitive (rank 3) to vertical hydraulic conductivity 
between layers 1 and 2 and between layers 4 and 5 


Mid-Hawthorn aquifer (layer 9) 


The aquifer is sensitive (rank 4) to: 


Vertical hydraulic conductivity between layers 5 and 6, layers 
6 and 7, layers 7 and 8, and layers 8 and 9 

Its own horizontal hydraulic conductivity 

Its own specific storage 

The boundary condition conductance 

River conductance 


Table 7-6. Significant sensitivity of the Water Table aquifer. 


Parameter Layer Multiplier Rank 
Kh 1 5 2 
Kh 1 0.2 2, 
Kh 1 0.1 2 
Kv 1 0.1 2 
Specific yield 1 0.5 2 
Specific yield 1 0.25 2 
Kh 3 10 2 
Kh 3 5 2 
Kh 5 10 2 
Kh 6 10 2 
Kh 6 5 2 
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Drain conductance N/A 0.1 2 
River conductance N/A 0.1 2 
River conductance N/A 0.01 2 
Evapotranspiration N/A 0.8 2 
Recharge N/A 1.2 2 
Recharge N/A 0.8 2 


Kh = horizontal hydraulic conductivity; Kv = vertical hydraulic conductivity; N/A = not applicable. 


Table 7-7. Significant sensitivity of the Lower Tamiami aquifer. 


Parameter Layer Multiplier Rank 
Kv 2 0.1 3 
Kh 3 10 3 
Kh 1 5 2 
Kh 1 0.2 2 
Kh 1 0.1 2 
Kv 1 0.2 2 
Kv 1 0.1 2 
Specific yield 1 0.25 2 
Kv 2 10 2 
Kv 2 5 2 
Kv 2 0.2 2 
Kh 3 5 2 
Kh 3 0.2 2 
Kh 3 0.1 2 
Drain conductance N/A 0.01 2 
River conductance N/A 100 2 
River conductance N/A 10 2 
River conductance N/A 0.01 2 
Evapotranspiration N/A 0.8 2 
Recharge N/A 1.2 2 
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Kh = horizontal hydraulic conductivity; Kv = vertical hydraulic conductivity; N/A = not applicable. 


Table 7-8. Significant sensitivity of the Sandstone aquifer — clastic zone. 


Parameter Layer Multiplier Rank 
Kv 2 0.2 4 
Kv 2 0.1 4 
Kv 4 0.2 4 
Kv 4 0.1 4 
Kh 5 0.1 4 
Kv 1 0.1 3 
Kv 4 10 3 
Kv 4 Э 3 
Kh 5 10 3 
Kh 5 0.2 3 
Kh 7 0.2 3 
Kh 7 0.1 3 
Kh 1 0.2 2 
Kh 1 0.1 2 
Kv 1 0.2 2 
Specific yield 1 0.5 2 
Specific yield 1 0.25 2 
Kh 3 10 2 
Kh 3 0.1 2 
Kh 5 5 2 
Kh 6 10 2 
Kh 6 2 2 
Ку 6 0.2 2 
Ку 6 0.1 2 
Kh 7 10 2 
Kh 7 5 2 
Drain conductance N/A 0.01 2 
River conductance N/A 10 2 
River conductance N/A 0.1 2 
River conductance N/A 0.01 2 
Evapotranspiration N/A 0.8 2 
Recharge N/A 1.2 2 


Kh = horizontal hydraulic conductivity; Kv = vertical hydraulic conductivity; N/A = not applicable. 
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Table 7-9. Significant sensitivity of the Sandstone aquifer — carbonate zone. 


Parameter Layer Multiplier Rank 
Kv 2 0.1 4 
Kv 4 10 4 
Kv 4 5 4 
Kv 4 0.2 4 
Kv 4 0.1 4 
Kv 1 0.1 3 
Kh 7 10 3 
Drain conductance N/A 0.01 3 
Kh 1 5 2 
Kh 1 0.2 2 
Kh 1 0.1 2 
Ку 1 0.2 2 
Kh 2 10 2 
Kv 2 0.2 2 
Kh 3 10 2 
Kh 3 0.2 2 
Kh 3 0.1 2 
Kh 5 10 2 
Kh 5 5 2 
Kh 5 0.1 2 
Kh 6 10 2 
Kh 6 5 2 
Kv 6 10 2 
Kv 6 5 2 
Ку 6 0.1 2 
Kh 7 5 2 
Specific storage 8 100 2 
Kh 9 10 2 
Kh 9 5 2 
Kh 9 0.1 2 
Drain conductance N/A 0.1 2 
River conductance N/A 100 2 
River conductance N/A 10 2 
River conductance N/A 0.1 2 
River conductance N/A 0.01 2 
Evapotranspiration N/A 0.8 2 
Recharge N/A 1.2 2 


Kh = horizontal hydraulic conductivity; Kv = vertical hydraulic conductivity; N/A = not applicable. 
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Table 7-10. Significant sensitivity of the Mid-Hawthorn aquifer. 


Parameter Layer Multiplier Rank 
Kv 5 0.2 4 
Kv 5 0.1 4 
Kv 6 0.2 4 
Kv 6 0.1 4 
Kv 7 10 4 
Kv 7 5 4 
Kv 7 0.2 4 
Kv 7 0.1 4 
Kv 8 10 4 
Kv 8 3 4 
Ky 8 0.2 4 
Ky 8 0.1 4 
Kh 9 10 4 
Kh 9 S 4 
Kh 9 0.2 4 
Kh 9 0.1 4 
Specific storage 9 100 4 
GHB conductance N/A 0.1 4 
GHB conductance N/A 0.01 4 
West GHB conductance N/A N/A 4 
River conductance N/A 0.1 4 
River conductance N/A 0.01 4 
Kv Б 10 3 
Ку 2 5 3 
Kv 6 10 3 
GHB conductance N/A 100 3 
Kh 1 5 2 
Kh 1 0.2 2 
Kh 1 0.1 2 
Specific yield 1 0.25 2 
Kh 2 10 2 
Kh 3 10 2 
Kh 4 10 2 
Kv 6 5 2 
Specific storage 9 0.1 2 
GHB conductance N/A 10 2 
North GHB conductance N/A N/A 2 
River conductance N/A 100 2 
River conductance N/A 10 2 
Evapotranspiration N/A 1.2 2 
Evapotranspiration N/A 0.8 2 
Recharge N/A 0.9 2 
Recharge N/A 0.8 2 


GHB = general head boundary; Kh = horizontal hydraulic conductivity; Kv = vertical hydraulic conductivity; N/A = not 
applicable. 
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7.4.1 Stresses and Related Variables 


The sensitivity scenarios included scaling back ET and recharge by factors of 0.8, 0.9, 1.1, and 1.2 versus 
the calibrated model. Four runs scaled ET and four runs scaled recharge. Table 7-11 presents the sensitivity 
results for the ET and recharge scenarios that exhibited at least slight sensitivity (rank 2). Of the ET 
scenarios, only 0.8- and 1.2-scaled runs exhibited slight model sensitivity, and only the 0.8-scaled run 
exhibited sensitivity in all aquifers. Extinction depth also was varied in four scenarios by the same scale 
factors; however, none of the scenarios exhibited even slight sensitivity (rank 2). Model sensitivity to gross 
recharge was slight (rank 2) for all aquifers. Gross recharge is the recharge input provided to MODFLOW 
and should not be confused with net recharge, which is the combined effect of gross recharge and ET. 


Table 7-11. Significant sensitivity to stresses. 


Aquifer Rank 
Parameter Multiplier WTA LTA SSA- SSA- MHA Rank Sum 
clastic | carbonate 
Evapotranspiration 0.8 2 2 2 2 2 10 
Evapotranspiration 1.2 1 1 1 1 2 6 
Recharge 1.2 2 2 2 2 1 9 
Recharge 0.9 1 1 1 1 2 6 
Recharge 0.8 2 1 1 1 2 7 


LTA = Lower Tamiami aquifer; MHA = Mid-Hawthorn aquifer; SSA = Sandstone aquifer; WTA = Water Table aquifer. 


7.4.2 Conclusion of the Sensitivity Analysis 


The model appeared reasonably calibrated; any change decreased model performance. The results from the 
global model sensitivity analysis showed that model-wide changes to any considered parameter did not 
improve the overall model calibration. Modifying the most sensitive parameters generally resulted in a net 
degradation of model performance. The sole exception (aside from perhaps miniscule improvements), 
discussed above (layer 6 Kh of 0.2x that of the calibrated model) results in the range of relative MAE 
extending upwards nearly 2 ft in the WTA. 


The calibration run was superior when considering the MAE across the model domain. While there may be 
instances where a sensitivity run might appear superior to the calibration run for one or more parameters at 
the local scale, or even for one or a few aquifers, this comes at the expense of poorer performance elsewhere. 


8 MODEL LIMITATIONS 


Overall, the transient groundwater flow model reasonably simulates regional hydrologic conditions across 
the model domain. Simulated potentiometric surfaces are in good visual agreement with regional 
potentiometric surface maps. Residual statistics also demonstrate good agreement between observed and 
simulated water levels. 


Although the LWCSIM adequately represents the groundwater flow patterns and fluxes within the model 
domain, there are a few areas where model residual head errors are significant, particularly in the Cape 
Coral area of model layer 9. Localized head changes near some major wellfields were difficult to simulate 
given the regional nature of the application. Furthermore, although the model generally reproduces temporal 
water level trends at most monitor wells, there are significant differences between observed and simulated 
temporal water level response at specific wells. These differences are not unexpected given the use of 
aquifer parameter zones during calibration and the limitations of available pumpage and hydrogeologic 
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data. In these instances, qualitative and quantitative calibration criteria were balanced to achieve the best 
overall calibration. Despite these differences, the calibrated model reasonably simulates the regional and 
three-dimensional hydrogeologic features of the model domain. 


Potential users of the LWCSIM should note because of recognized data limitations, model simulation is 
more appropriate at sub-regional and regional scales rather than very local or site-specific scales for 
simulation of hydrologic conditions. Grid cell sizes are a uniform 1,000 ft. Any simulation results at a scale 
of a few hundred feet would be problematic unless a finer grid discretization is used. The conceptual model 
used to construct the model is a simplified representation of the true groundwater flow system. The SAS 
and IAS in the LWC Planning Area can be characterized as a complex, heterogeneous aquifer system. 


Various model parameters and boundary conditions were determined based on land use distribution, 
configuration of waterways and other surface water bodies, and climatic conditions during the calibration 
period. If any of these factors were to significantly change, parameters such as recharge, ET, and river 
conductance may need to be modified. 


9 RECOMMENDATIONS FOR FUTURE DATA COLLECTION 
9.1 Water Use Data 


Currently, the SFWMD has a good data reporting system in place for PS water use. For a variety of reasons, 
the reported water use data for AG, L/R, and DSS are less reliable. Therefore, modelers use various tools 
to estimate water use from each category and eventually use those as model calibration parameters. Better 
data for these water use categories would result in greater confidence in model calibration to support water 
supply planning efforts. 


9.2 Water Level (Head) Data 


The assessment of model performance is based on available data. At many stations where data are available, 
there are months with no data points, which makes it difficult, and sometimes impossible, to assess model 
performance at a local scale or over short time periods (e.g., a single dry season). It probably is more 
valuable to have a moderately accurate measure of mean head over a time interval than to have a highly 
accurate measure at one instant. One recommendation is to collect groundwater data at more locations, 
especially in places that are well connected to withdrawal areas and are far away from areas where data 
already are collected. In addition, it may be more valuable to collect more groundwater data during times 
of extreme water levels and less during the transitions between extremes. Again, better data to support 
model calibration, especially during extreme events, result in greater confidence in the model's predictive 
capabilities. 


9.3 Hydraulic Conductivity Data 


In general, there is a good coverage of field data for all the aquifers to determine hydraulic properties, 
especially where the aquifers are heavily used. However, where the aquifers are not heavily used, there is 
a scarcity of field data, especially in the SSA and MHA. It is essential to characterize system hydraulic 
properties as accurately as possible when building a regional groundwater model. Due to the complex 
heterogeneity of the subsurface system, interpolating hydraulic properties from known areas to areas where 
the field measurements are not available can introduce significant errors in the model. Therefore, it is 
recommended that additional wells be installed and aquifer tests be conducted in the SSA and MHA of the 
LWC Planning Area. 
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APPENDIX А: 
NEXRAD CORRECTION 


INTRODUCTION 


Rainfall plays a prominent role in the development and application of hydrologic models and water budgets. 
Traditionally, rain gauges have been used to measure and estimate the rainfall input needed for hydrologic 
models and water budgets. However, errors associated with rain gauge measurements limit their use and 
increase the uncertainties in the results (Habib et al. 2009, Skinner et al. 2009, Craciun and Catrina 2016). 
In addition, estimating rainfall over large areas requires maintaining a large network of rain gauge stations, 
which can be costly and limited by accessibility to remote station locations. Advances in radar technology 
and continuous improvement in rainfall volume detection and estimation have made radar-derived rainfall 
data useful for hydrologic applications (Habib et al. 2007, 2009; Skinner et al. 2009; Craciun and Catrina 
2016). 


In the early 1990s, the National Weather Service established the Next-Generation Weather Radar 
(NEXRAD) system (approximately 160 ground-based radar stations) across the United States, including 
five stations in Florida (Habib et al. 2009, Skinner et al. 2009). Unlike rain gauge networks, NEXRAD can 
estimate rainfall over a large area by providing a better representation of rainfall distribution (Skinner et al. 
2009). However, there are multiple forms of error associated with radar-derived data, including systematic 
bias, random error, and temporal and spatial radar error (Habib et al. 2007). These, errors greatly affect the 
accuracy of rainfall estimates; therefore, in most cases, evaluation and correction of radar data are needed 
to reduce uncertainties. Comparisons of rain gauge and radar estimates indicate radar tends to overestimate 
small rain volumes and underestimate large rain volumes (Over et al. 2007, Hardegree et al. 2008, Skinner 
et al. 2009). 


Improvement of radar-derived data can be achieved by correcting for inherent bias using adjustment factors; 
however, systematic deviations in radar data from reference surface rainfall (rain gauges) can remain after 
correction (Habib et al. 2007). Several approaches for correcting NEXRAD data have been investigated 
and implemented successfully depending on period, resolution, and area of interest (Chumchean et al. 2006, 
Hanchoowong et al. 2012, Vernimmen et al. 2012, Smith and Rodriguez 2017). 


The South Florida Water Management District (SFWMD) uses NEXRAD products in the development of 
water budgets and application of hydrologic models. However, there are concerns about the accuracy of the 
NEXRAD rainfall data, which is a key component in these model applications (Beeson et al. 2011). 
Multiple evaluations have been completed to address these concerns (Skinner et al. 2009, Brown 2014, 
Trimble 2018). 


The objectives of this study were to 1) evaluate the relationship between rain gauge and NEXRAD rainfall 
estimates during a 16-year period (1999 to 2014) in the Lower West Coast (LWC) Planning Area of South 
Florida; 2) determine the amount of error in the NEXRAD data set; and 3) develop a simple systematic 
approach for correcting daily NEXRAD rainfall using adjustment factors derived from the bias between 
rain gauge and NEXRAD rainfall. In this study, rain gauges were assumed to offer a better representation 
of rainfall accumulation than NEXRAD at specific locations because they are direct measurements. 
Therefore, rain gauges were used as a ground truth for rainfall despite the presumed inaccuracies and 
limitations of the rain gauge data set. 
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METHODS 
Study Area 


The LWC Planning Area encompasses approximately 5,130 square miles in the southwestern region of 
Florida between the Gulf of Mexico and the western edge of the Everglades, extending north to the 
Southwest Florida Water Management District jurisdictional boundary and south to include portions of 
Everglades National Park (Figure A-1). The LWC Planning Area average elevation is approximately 
16 feet National Geodetic Vertical Datum of 1929 (NGVD29), with elevations ranging from slightly above 
sea level to as high as 65 feet NGVD29 (SFWMD 2017). Most rainfall in the LWC Planning Area occurs 
during the wet season, which normally is from June through September; the remaining months typically are 
considered dry months. The amount of rain varies on a monthly and yearly basis and is greatly affected by 
climatic disturbances such as regional droughts, seasonal storm events, and hurricanes. In addition, rainfall 
can vary spatially depending on the geographic features of the area. Rainfall estimates in the LWC Planning 
Area were derived using the extensive rain gauge network maintained by the SFWMD (Figure A-1) and/or 
using NEXRAD. 


Rain Gauge and NEXRAD Data 


Rainfall data in the LWC Planning Area are collected by a network of rain gauges maintained by the 
SFWMD. The number of active rain gauges varied from 65 to 82 gauges during the study period (1999 to 
2014). During the study period, the SFWMD operated four rain gauge recording types: tipping buckets, 
weighing buckets, float-type, and simple standard. Long-term variability in the number of available rain 
gauges was due to optimization of the monitoring network, while short-term changes likely were associated 
with equipment maintenance and/or malfunction. Therefore, missing data values are common over time. A 
comprehensive quality assurance and quality control (QA/QC) process was conducted by the SFWMD to 
address possible data problems before data were loaded into the agency's DBHYDRO database. In addition, 
rain gauge data quality and accuracy can be affected by wind speed, canopy coverage, and physical 
obstructions that may be difficult to discern during the QA/QC process. However, a 2014 QA/QC 
evaluation indicated rain gauge measurement data better represents rainfall amounts at co-located rain 
gauge and NEXRAD pixels (Brown 2014). Daily rain gauge values generally are the result of daily sum 
values; therefore, finer time resolution may be available at some rain gauge stations 


The SFWMD acquires NEXRAD data as part of an agreement with a private vendor, who provides data to 
all five water management districts within the State of Florida. NEXRAD data are obtained in near real 
time at 15-minute intervals, which is aggregated into hourly and daily values. Prior to delivery, the 
near-real-time data are adjusted by the vendor using selected rain gauge data, which are incorporated into 
the algorithm computation (Brown 2014) and used to validate and calibrate real-time radar rainfall data 
(Craciun and Catrina 2016). NEXRAD data are delivered to the SFWMD in a 2 x 2-kilometer grid format 
that includes a 35-mile buffer around the SFWMD’s boundary; however, database queries only provide 
access to pixel data within the SFWMD’s jurisdictional boundary. Each 2 x 2-kilometer grid is equivalent 
to 1 pixel, which is identified by a unique pixel identification number in DBHYDRO (Figure А-2). For this 
study, a grid of 4,763 pixels covering the entire LWC Planning Area was used. All rain gauge and NEXRAD 
data used in this evaluation were obtained from DBHYDRO. 
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Figure A-1. Rain gauge stations within the Lower West Coast Planning Area. 
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Figure A-2. Florida radar locations and the Next-Generation Weather Radar (NEXRAD) rain grid along 
with corresponding pixels in reference to rain gauge stations in the Lower West Coast 
Planning Area. 
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Statistical and Spatial Methods 


Each rain gauge and associated 2 x 2-kilometer pixel were compared on a daily, monthly, and yearly basis 
from 1999 to 2014. Mean rainfall volume difference for monthly and yearly totals was calculated as follows: 


RD = NRD; — RG; (1) 


Where RD is the mean rainfall difference, and NRD, and RG; represent the mean rainfall estimated by 
NEXRAD and recorded by the rain gauge, respectively. The relative bias between the monthly and yearly 
recorded rain gauge and estimated NEXRAD rainfall totals was calculated as follows: 


Relative bias = — x 100 (2) 
i 


Linear regression of yearly, monthly, and daily rain gauge and NEXRAD rainfall totals was performed to 
evaluate the overall relationship between the variables. Furthermore, the Kolmogorov-Smirnov statistical 
test for normality showed a marked positive skew of the daily values, indicating a considerable number of 
small values near or equal to zero; this observation was not surprising because it does not rain every day. 
Therefore, the relationship between NEXRAD data and individual rain gauges was evaluated using 
Kendall’s tau (t) rank-based correlation method (Helsel and Hirsch 1992). Monthly correlation for each 
rain gauge and NEXRAD pixel pair was calculated using daily values (n = 14), and correlation coefficients 
were used to discern which rain gauges displayed a good relationship (t > 0.75) with the NEXRAD pixels. 
Rain gauges showing a poor relationship with the NEXRAD pixels were excluded from further analysis. 
Yearly and monthly rain gauge and NEXRAD totals were calculated. The bias between the rain gauge and 
NEXRAD was determined by the ratio between the two (Steiner et al. 1999, Goudenhoofdt and Delobbe 
2009, Smith and Rodriguez 2017). The resulting value was used as the NEXRAD adjustment factor, 
calculated as follows: 


AdjF (bias) = ул /у + 3) 


Where AdjF is the multiplicative factor used to adjust the daily NEXRAD value, and RG; and МЕР); are the 
monthly rainfall recorded by the rain gauge and estimated by NEXRAD, respectively, during an equally 
match-paired number of days (n > 14). To prevent overcorrection of NEXRAD daily values and to include 
most adjustment factors, the adjustment factor was arbitrarily capped at not greater than three. Therefore, 
daily NEXRAD values with adjustments factors greater than three received a maximum adjustment factor 
of three. Visual examination of the resulting bias showed several rain gauges located less than 2 kilometers 
from each other, which can produce inconsistencies in the bias arising from high and low values in spatially 
close stations. Therefore, the resulting monthly bias value at each rain gauge station was normalized using 
а complete hierarchical cluster analysis (Müllner 2013, Máchler et al. 2018), in which bias points were the 
result of the averaged bias of each rain gauge located within 2 kilometers of each other. The resulting 
monthly bias was used as the multiplicative adjustment factor to correct the daily NEXRAD pixel value 
using geostatistical processes. 


Adjustment factors were interpolated across the NEXRAD grid using the automated autoKrige method 
described by Hiemstra et al. (2009). The autoKrige function tests different interpolation models (linear, 
spherical, exponential, Gaussian, and Matern) by estimating semi-variograms and selecting the best fit 
kriging model (Hiemstra et al. 2009). This process produced a monthly bias grid with adjustment factors at 
each NEXRAD pixel; these adjustment factors were applied to the uncorrected NEXRAD values to produce 
corrected NEXRAD values. Point values for the bias and the uncorrected and corrected NEXRAD were 
converted into a raster for visual evaluation. During this process, a few inconsistencies associated with the 
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NEXRAD and rain gauge pairing were observed in the corrected МЕХКАР raster. These inconsistencies 
resulted from mismatched NEXRAD and rain gauge daily values. The data were examined and removed 
from daily values. Subsequently, the entire process described above was repeated and a new set of raster 
plots was created. The difference or change in rain was calculated by subtracting the corrected from the 
uncorrected NEXRAD, and the results were displayed in a separate raster. All statistical and spatial data 
analyses were performed using R programming language (version 3.5.1; R Development Core Team 2018). 


Previous analysis of the NEXRAD data set reported missing rainfall for April 2001 (Brown 2014). To 
maintain a continuous record, the missing rainfall data were produced using the methods described in the 
documentation for the South Florida Water Management Model version 5.5 (SFWMD 2005). 


RESULTS 
Comparison of Rain Gauge and Uncorrected NEXRAD 


Rain gauge and estimated NEXRAD monthly average rainfall for the LWC Planning Area from 1999 to 
2014 followed similar trends, with slight differences in magnitude associated with periods of higher rainfall 
volumes (Figure A-3a). Yearly comparisons of monthly averages for rain gauge and estimated NEXRAD 
rainfall showed marked differences in values reported from 2002 through 2004 (Figure A-3b). The mean 
rainfall volume difference between the rain gauge and NEXRAD monthly totals showed a large 
underestimation of rainfall by NEXRAD from 2002 to 2004; this trend was less noticeable but consistent 
through 2007. While 2008 and 2009 showed a consistent overestimation of NEXRAD monthly rainfall 
totals (Figure A-4), subsequent years showed greater agreement between the rain gauge and NEXRAD 
monthly totals. 


Scatter plots comparing daily, monthly, and yearly rainfall totals between the rain gauges and associated 
uncorrected NEXRAD pixels indicated a better fit or higher relationship of the data using monthly totals 
(Figure A-5). The monthly coefficient of determination (denoted by R?) obtained from the linear regression 
for the entire period was greater than the yearly and daily coefficients. Further comparison of regression 
coefficients calculated from the daily, monthly, and yearly totals confirmed that monthly relationships were 
considerably better than yearly and daily relationships for all years. However, the daily relationships 
between rain gauge and NEXRAD data improved from 2007 through 2014 (Table A-1). Monthly totals 
were determined to be more appropriate in the calculation of bias between the rain gauge and corresponding 
NEXRAD pixel. In addition, monthly bias was expected to be suitable to capture seasonal changes and 
climate disturbances such as droughts and storm events through the years. To ensure appropriate adjustment 
of the daily rainfall NEXRAD values, the calculated bias was restricted using correlation coefficients 
(т = 0.75) derived from the monthly rain gauge and NEXRAD relationship. The restriction was based on 
the potential for using inaccurately paired rain gauge and NEXRAD values to calculate the bias due to 
inherent errors that exist in both methods. Therefore, after the correlation restriction, the number of rain 
gauge stations and associated bias used for the NEXRAD correction was reduced by 3246 from 1999 to 
2014. Excluding these rain gauge stations during the spatial interpolation reduced the potential transfer of 
errors across the entire area. In addition, higher correlations on paired rain gauge and NEXRAD values may 
be the result of prior adjustment to the NEXRAD data by the vendor when incorporating rain gauge data 
into the algorithm calculation. 
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Figure A-5. Scatter plots comparing all paired rain gauge and estimated uncorrected Next-Generation 
Weather Radar (NEXRAD) for (a) yearly, (b) monthly, and (c) daily rainfall totals, in 
inches, from 1999 to 2014. 


Table А-1. Yearly coefficients of determination (R°) calculated from the relationship between all 
paired rain gauges and uncorrected Next-Generation Weather Radar (NEXRAD) data using 
yearly, monthly, and daily totals. The number of observations (n) indicates the number of 


rain gauge and uncorrected NEXRAD pairs available for each year. 
Year Yearly Totals Monthly Totals Daily Totals 
n R? n R? n R? 

1999 75 0.63 740 0.82 19,461 0.49 
2000 74 0.79 690 0.85 19,496 0.6 
2001 76 0.79 700 0.88 19,753 0.5 
2002 73 0.42 748 0.8 19,935 0.53 
2003 79 0.81 797 0.85 21,576 0.57 
2004 79 0.7 813 0.84 23,206 0.57 
2005 77 0.61 856 0.81 22,685 0.56 
2006 76 0.73 783 0.86 22,801 0.59 
2007 82 0.78 841 0.84 23,222 0.6 
2008 81 0.67 889 0.91 25,579 0.76 
2009 79 0.67 882 0.91 25,544 0.71 
2010 81 0.65 928 0.89 26,133 0.7 
2011 81 0.73 859 0.9 24,894 0.69 
2012 65 0.7 726 0.9 21,513 0.74 
2013 65 0.6 742 0.92 21,432 0.74 
2014 66 0.57 756 0.88 21,106 0.73 
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Uncorrected versus Corrected NEXRAD 


Comparison of monthly totals for selected paired rain gauge and uncorrected and bias-corrected NEXRAD 
data between 1999 and 2014 showed an improvement in the coefficient of determination (R2) and 
relationship (Figure A-6). Yearly regression coefficients calculated from monthly totals showed similar 
improvement after bias correction was applied to the daily NEXRAD values (Table A-2). In addition, 
yearly rain gauge and corrected NEXRAD rainfall totals were more equivalent than uncorrected NEXRAD 
totals (Figure А-7). The mean yearly relative bias of monthly totals between the rain gauge and uncorrected 
NEXRAD data indicated that, on average, NEXRAD underestimated rainfall from 1999 to 2003 by as much 
as 2396. Subsequent years showed slight (1% to 7%) overestimation and underestimation of NEXRAD 
yearly totals (Figure A-8). The percent error between the rain gauge and NEXRAD data was greatly 
reduced (less than 596) for all years after the correction was applied to the daily NEXRAD values. When 
NEXRAD underestimated rainfall, the corrected NEXRAD data resulted in a slight overestimation 
compared to rain gauge yearly values, which suggests a possible overcorrection of the NEXRAD daily 
values (Figure A-8). 
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Figure A-6. Scatter plots comparing selected paired rain gauge and estimated Next-Generation Weather 
Radar (NEXRAD) monthly rainfall totals, in inches, for daily uncorrected and corrected 
NEXRAD values from 1999 to 2014. Circles represent a monthly value for each rain gauge 
and NEXRAD pixel located within the Lower West Coast Planning Area. 
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Table A-2. Yearly coefficients of determination (R°) calculated from the relationship between selected 
paired rain gauges and Next-Generation Weather Radar (NEXRAD) monthly rainfall totals 
before and after correction. The number of observations (n) indicates the number of rain 
gauge and NEXRAD monthly pairs available for each year with a good correlation 


(т > 0.75). 
Uncorrected NEXRAD Bias-Corrected NEXRAD 
Year 
n R? n В? 
1999 740 0.82 752 0.85 
2000 690 0.85 696 0.88 
2001 700 0.88 701 0.88 
2002 748 0.8 755 0.86 
2003 797 0.85 807 0.87 
2004 813 0.84 816 0.88 
2005 856 0.81 859 0.86 
2006 783 0.86 787 0.9 
2007 841 0.84 845 0.89 
2008 889 0.91 893 0.93 
2009 882 0.91 880 0.94 
2010 928 0.89 927 0.96 
2011 859 0.9 867 0.91 
2012 726 0.9 725 0.93 
2013 742 0.92 745 0.94 
2014 756 0.88 756 0.9 
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Figure А-7. Comparison of paired rain gauge and estimated Next-Generation Weather Radar 
(NEXRAD) yearly rainfall totals, in inches, for uncorrected and corrected NEXRAD values 
from 1999 to 2014. 
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Figure A-8. Mean yearly relative bias of monthly totals between available rain gauges and the 
uncorrected and corrected Next-Generation Weather Radar (NEXRAD). Error bars 
represent + 1 standard error of the mean. 


Bias 


The monthly bias, or adjustment factor, is defined by the sum of daily rainfall values recorded by the rain 
gauge divided by the sum of the associated estimated daily rainfall NEXRAD value. Monthly bias values 
greater than 1 indicated underestimation, while bias values lower than 1 suggested overestimation of rainfall 
by NEXRAD. Monthly bias derived from the difference between the rain gauge and NEXRAD rainfall 
resulted in approximately 2.7% of the bias greater than 3, suggesting a gross underestimation of NEXRAD 
rainfall. When the bias value was greater than 3, a maximum bias or adjustment factor of 3 was applied to 
the daily NEXRAD value to prevent overcorrection. Bias values for the uncorrected NEXRAD appeared 
normally distributed, suggesting a tendency to equally overestimate and underestimate monthly rainfall 
totals (Figure A-9a). As bias is calculated monthly and applied to the daily NEXRAD value, some bias is 
anticipated to remain in the corrected NEXRAD monthly totals. Evaluation of the remaining monthly bias 
after NEXRAD rainfall correction showed a substantial reduction in the differences between the rain gauge 
and NEXRAD data (Figure A-9b). Overall, the 1999 to 2014 monthly bias after correction of the NEXRAD 
rainfall showed a statistically tighter proximity to 1 compared to the uncorrected NEXRAD bias 
(Table A-3). Further evaluation of the bias showed an improvement of the uncorrected NEXRAD data after 
2007. This trend also was highlighted in the previous comparison of uncorrected and corrected NEXRAD 
rainfall yearly totals. Finally, continuous improvement of the NEXRAD correction methodology by the 
vendor has resulted in better estimation of NEXRAD rainfall data compared to previous years. 
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Figure A-9. Density plot comparison of monthly bias from 1999 to 2014 at each rain gauge station 
(a) before and (b) after Next-Generation Weather Radar (NEXRAD) rainfall correction. 
Values greater or lower than 1 represent overestimation or underestimation, respectively, of 


NEXRAD rainfall monthly values, in inches. 


Table A-3. Summary statistic comparison of monthly bias from 1999 to 2014 for uncorrected and 
corrected Next-Generation Weather Radar (NEXRAD) rainfall. 


Же; В1а$ 
о Uncorrected NEXRAD Rainfall Corrected NEXRAD Rainfall 
Maximum 3.000 2.978 
90" percentile 1.667 1.351 
75 percentile 1.224 1.079 
Mean 1.095 1.007 
Standard Deviation 0.529 0.354 
Standard Error Mean 0.005 0.003 
Median 1.000 0.994 
25% percentile 0.811 0.853 
10" percentile 0.596 0.644 
Minimum 0.006 0.007 


Spatial Analysis 


The monthly bias, or adjustment factors, derived from the relationship between the rain gauge and 
NEXRAD rainfall were interpolated across the LWC Planning Area to produce 192 raster plots, one for 
each month between 1999 and 2014. The interpolated NEXRAD bias did not show any spatial or temporal 
trends in the distribution and variability of the bias. Instead, the variability of the bias appeared to be 
influenced primarily by the relationship between the rain gauge and NEXRAD pixel. Results showed large 
variability in the number of rain gauges available for each month, which was heavily influenced by the rain 
gauge data quality, the amount of rainfall recorded, and the number of paired rain gauge and NEXRAD 
observations available for each month. In addition, the southernmost portion of the LWC Planning Area 
has a lack of rain gauge monitoring stations, whereas the middle to northern region has a larger 
concentration and clustering of rain gauge monitoring stations. Overall, the NEXRAD bias appeared 
localized in relation to rain gauge locations. 


As described earlier, the clustering of rain gauge stations was minimized by normalizing adjustment factors 
that appeared within 2 kilometers of each other. The remaining adjustment factors were used to correct the 
daily NEXRAD pixel values (Figure A-10). Two additional raster plots were produced to visually compare 
the uncorrected and corrected NEXRAD monthly totals. In some cases, differences between the uncorrected 
and corrected NEXRAD values were visually discernable, particularly during months with higher rainfall 
(Figure A-10a). However, months with low rainfall did not appear to visually differ (Figure A-10b). Large 
variability in monthly rainfall volumes occurred across the LWC Planning Area, suggesting heavy 
downpours in localized areas. The change in monthly rainfall total was spatially depicted by the difference 
between the uncorrected and corrected NEXRAD values. The resulting raster plot (denoted by Change in 
NEXRAD) showed where the monthly NEXRAD total rainfall in inches was increased or reduced above 
the zero mark (Figure A-10). 


As shown in the monthly total raster plots, one advantage of applying monthly bias adjustments to the daily 
NEXRAD values is the process can capture seasonal changes and climate disturbances such as droughts 
and storm events through the years. This is particularly important for periods associated with tropical storms 
and hurricanes when rain gauge data may be limited or nonexistent in areas severely affected by the storm. 
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Figure A-10. Lower West Coast Planning Area raster maps showing the calculated monthly bias, monthly rainfall totals (in inches) for 
uncorrected and corrected Next-Generation Weather Radar (NEXRAD), and net change in NEXRAD after the adjustment is applied 
for (a) a dry season month in December 2003 (X200312), and (b) a rainy season month in July 2006 (X200606). 
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SUMMARY AND CONCLUSION 


The main objectives of this study were to 1) evaluate the relationship between rain gauge measurements 
and estimated NEXRAD values using yearly, monthly, and daily rainfall totals in the LWC Planning Area 
of Florida during a 16-year period (1999 to 2014); and 2) develop a method to correct NEXRAD rainfall 
bias using adjustment factors derived from the relationship of co-located rain gauge and NEXRAD pixel 
pairs. An important assumption in this study was that rain gauge measurements offer a better representation 
of rainfall accumulation than estimated NEXRAD values at the same location. Although, the quality and 
continuity of the rain gauge measurement data can be affected by multiple factors, previous QA/QC 
evaluations indicated rain gauge measurements are a better representation of rainfall than NEXRAD values 
at co-located rain gauge and NEXRAD pixels (Brown 2014). 


Overall, NEXRAD underestimated rainfall in the LWC Planning Area from 1999 to 2003 by as much as 
23%, while only slight (1% to 7%) overestimations and underestimations were observed in subsequent 
years. Results showed improvement in the quality of NEXRAD data provided by the vendor from 2007 
through 2014, as indicated by higher coefficients of determination derived from the rain gauge and 
NEXRAD daily and monthly rainfall totals. 


Comparison of rain gauge measurements and estimated NEXRAD values showed a stronger correlation 
using monthly totals than yearly and daily totals. As a result, monthly correction factors were determined 
to be better suited for the correction of daily NEXRAD values. Furthermore, adjustment of daily NEXRAD 
values using monthly bias factors provided a way to capture seasonal changes and climate disturbances 
such as droughts and storm events. This was particularly important for periods associated with tropical 
storms and hurricanes when rain gauge data appeared affected, limited, or nonexistent in areas severely 
affected by the storm. 


Correction of the NEXRAD daily rainfall data was achieved using a systematic approach in which 
calculated monthly adjustment factors were interpolated across the LWC Planning Area and applied to the 
daily values. Although corrections of daily NEXRAD data were applied systematically to all 4,763 pixel 
values in the LWC Planning Area, there were minor inconsistencies in daily rainfall values at co-located 
rain gauge and NEXRAD pixel pairs that could have resulted in localized overestimation or underestimation 
of rainfall. Therefore, further examination of the raster plots was necessary to find and manually remove 
inconsistencies before producing the final adjusted NEXRAD data set. Rather than completely eliminating 
the bias between the rain gauge and NEXRAD values, the process developed in this study attempted to 
reduce the bias while errors associated with rain gauges were preserved. 


After applying the monthly correction to the daily NEXRAD value, differences in rain gauge and NEXRAD 
daily, monthly, and yearly rainfall totals were reduced. The monthly relationship of co-located rain gauge 
and NEXRAD pixel pairs increased after the correction, while monthly bias between the rain gauge and 
NEXRAD decreased. 


Finally, results showed rain gauge data can be used to improve the accuracy of estimated NEXRAD daily 
rainfall volumes. In addition, the process for correcting large NEXRAD data sets was optimized and 
partially automated using R scripts. This process reduced the time and effort needed to evaluate the 
relationship between the rain gauge and NEXRAD rainfall, allowed for the use of yearly and monthly 
adjustment factors, and expedited the production of corrected NEXRAD data and raster plots. The method 
described in this study can be applied to different areas of the SFWMD. The corrected NEXRAD data can 
be used as an alternative data set for the application and development of hydrologic models and water 
budgets. 
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APPENDIX В: 
LAYER TOPS/BOTTOMS AND THICKNESSES 
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Figure B-1. Topographic surface elevation of the Water Table aquifer. 
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Figure B-2. Bottom elevation of the Water Table aquifer. 
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Figure B-3. Bottom elevation of the Tamiami confining unit. 
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Figure B-4. Bottom elevation of the Lower Tamiami aquifer. 
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Figure B-5. Bottom elevation of the Upper Hawthorn confining unit. 
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Figure B-6. Bottom elevation of the Sandstone aquifer — clastic zone. 
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Figure B-7. Bottom elevation of the Sandstone aquifer confining unit. 
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Figure B-8. Bottom elevation of the Sandstone aquifer — carbonate zone. 
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Figure B-9. Bottom elevation of the Mid-Hawthorn confining unit. 
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Figure B-10. Bottom elevation of the Mid-Hawthorn aquifer. 
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Figure B-12. Thickness of the Tamiami confining unit. 
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Figure B-13. Thickness of the Lower Tamiami aquifer. 
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Figure B-14. Thickness of the Upper Hawthorn confining unit. 
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Figure B-15. Thickness of the Sandstone aquifer - clastic zone. 
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Figure B-16. Thickness of the Sandstone aquifer confining unit. 
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Figure B-17. Thickness of the Sandstone aquifer — carbonate zone. 
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Figure B-18. Thickness of the Mid-Hawthorn confining unit. 


B-19 


DESOTO 


SARASOTA | | ; a 


LAKE 
OKEECHOBEE 


Charlotte #8 
Harbor @ 


HENDRY 
COLLIER 


ROWARD 


———— 


Layer Thickness Original 
Layer 9 (Mid-Hawthorn Aquifer) 
Feet 
mm High : 268 


== Low:0 


ЕЭ LWwCSIM Model Extent 20 Miles 


А 4 г... 
C LWCSIM Active Model Boundary 2 20 Klometere 
= Ëi 3 ЖАМ 
\\ad.sfwmd.gov\dfsroot\G S\GSBiz\WS\LWCSIM\Report_Figures\mxd\ThicknessOrig_Lyr9.mxd 


Figure B-19. Thickness of the Mid-Hawthorn aquifer. 
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Figure B-20. Simulated bottom elevation of the Water Table aquifer. 
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Figure B-21. Simulated bottom elevation of the Tamiami confining unit. 
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Figure B-22. Simulated bottom elevation of the Lower Tamiami aquifer. 
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Figure B-23. Simulated bottom elevation of the Upper Hawthorn confining unit. 
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Figure B-24. Simulated bottom elevation of the Sandstone aquifer — clastic zone. 
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Figure B-25. Simulated bottom elevation of the Sandstone aquifer confining unit. 
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Figure B-26. Simulated bottom elevation of the Sandstone aquifer — carbonate zone. 
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Figure B-27. Simulated bottom elevation of the Mid-Hawthorn confining unit. 
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Figure B-28. Simulated bottom elevation of the Mid-Hawthorn aquifer. 
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Figure B-29. Simulated thickness of the Water Table aquifer. 
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Figure B-30. Simulated thickness of the Tamiami confining unit. 
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Figure B-31. Simulated thickness of the Lower Tamiami aquifer. 
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Figure B-32. Simulated thickness of the Upper Hawthorn confining unit. 
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Figure B-33. Simulated thickness of the Sandstone aquifer — clastic zone. 
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Figure B-34. Simulated thickness of the Sandstone aquifer confining unit. 
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Figure B-35. Simulated thickness of the Sandstone aquifer — carbonate zone. 
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Figure B-36. Simulated thickness of the Mid-Hawthorn confining unit. 


B-37 


DESOTO 


SARASOTA 


LAKE 
OKEECHOBEE 


Charlotte 
Harbor 


HENDRY 


~ “COLLIER 


TT BROWARD 


Layer Thickness 30-ft 
Minimum Layer Thickness 


Layer 9 (Mid-Hawthorn 
Aquifer) 

Feet 

шт High : 268 


== Low: зо 


ЕЭ LWCSIM Model Extent = | 20 мез 
L—1 LWCSIM Active Model Boundary р 


30 Kilometers 


Аз MAN 
\\ad.sfwmd.gov\dfsroot\GS\GSBiz\WS\LWCSIM\Report_Figures\mxd\Thickness30ftmin_Lyr9.mxd 


Figure B-37. Simulated thickness of the Mid-Hawthorn aquifer. 
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APPENDIX С: 
DOMESTIC SELF-SUPPLY COVERAGE PROCEDURE 


An important part of the Lower West Coast Surficial and Intermediate Aquifer Systems Model (LWCSIM) 
is the placement of pumping wells via the WEL package. The WEL package includes wells representing 
the full range of water being extracted from the various shallow and intermediate aquifers throughout the 
Lower West Coast (LWC) Planning Area. Public supply wells, agricultural and commercial irrigation wells, 
and domestic self-supply (DSS) wells are represented in the model. This memo describes the compilation 
of DSS wells for the WEL package. 


DSS well data from Lee, Collier, Glades, and Hendry counties were compiled from several sources. The 
initial source was the Marco Water Engineering, Inc. (2006) groundwater flow model of the LWC Planning 
Area surficial aquifer system. Chapter 2 of the Marco Water Engineering, Inc. (2006) report explained that, 
while residential water use represents the smallest demand compared to other categories, the projected 
population increase in coastal regions (primarily in Lee and Collier counties) could significantly impact 
(increase) groundwater demand. Due to limited funding and well data availability, Marco Water 
Engineering, Inc. (2006) only considered the 18,432 DSS wells in Lee County (Figure C-1). 


The Marco Water Engineering, Inc. (2006) model report did not detail how the DSS well coverage was 
derived. Therefore, the South Florida Water Management District (SFWMD) chose to compile the 
LWCSIM DSS well coverage using a combination of internal (SFWMD) and external (county) information 
sources. The SFWMD requested information on DSS wells, including location, casing depth, and total 
depth, from county officials in Lee, Collier, Hendry, and Glades counties. DSS well information was 
provided as geographic information system (GIS) shapefiles or Excel spreadsheets that were converted to 
GIS shapefiles so all DSS wells could be located in the model domain. Information from Lee and Collier 
counties was the most extensive and most complete. County DSS well data were combined into a 
comprehensive, consistent Excel table containing all pertinent data fields, including county, state plane 
coordinates, well casing depth, well total depth, and per capita use rate (PCUR). GIS was used to intersect 
the well locations with the LWCSIM grid to assign a row and column to each well. 


The PCUR initially was calculated from information in the 2012 LWC Water Supply Plan Update (SFWMD 
2012). The plan update contains projections of DSS populations and water demand by county. Specifically, 
in Appendix A of the plan update, Tables A-1 to A-8 contain the data that were used to calculate a PCUR 
for all DSS wells in Lee, Collier, Hendry, and Glades counties. According to the plan update, the PCUR is 
the total annual water use divided by the permanent residents. For example, to calculate the net average 
PCUR for DSS wells in Lee County in 2005: 


1. DSS population = 68,566 permanent residents 
2. Net average DSS water demand = 8.37 million gallons per day (mgd) = 8,370,000 gallons per day 


(gpd) 
3. PCUR = 8,370,000 gpd + 68,566 permanent residents = 122 gpd per person (gpd/p) 


Тһе PCUR also was calculated for net 1-in-10-year drought conditions (129 gpd/p), gross average 
(147 gpd/p), and gross 1-in-10-year drought conditions (156 gpd/p). The average of those four PCURs is 
139 gpd/p. 


Each county's overall average PCUR based on the 2012 LWC Water Supply Plan Update are as follows: 


Collier County = 219 gpd/p 
Glades County = 139 gpd/p 
Hendry County = 131 gpd/p 
Lee County = 139 gpd/p 
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Figure C-1. Domestic Self-Supply well coverage (From: Marco Water Engineering, Inc. 2006). 
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Discussions with SFWMD Water Supply Planning staff indicated a better estimate of the PCUR might be 
obtained using the SFWMD’s Water Supply Utility Project (WaSUP) database of finished residential water 
demand as reported by each public supply utility. Based on these numbers, each county’s average PCUR 
was adjusted to the following: 


Collier County = 123.67 gpd/p 
Glades County = 114.58 gpd/p 
Hendry County = 67.43 gpd/p 
Lee County = 71.05 gpd/p 


The derived PCURs were used as starting points for assigning extraction rates to model cells representing 
DSS wells. Rates per model layer per model cell were determined and assigned using well casing depth and 
total depth to determine which model layer(s) were being tapped, accounting for model cells that contain 
more than one DSS well. It was decided that the PCUR represents water used per person and that each 
household well provides water for an average of three persons. Therefore, a factor of three was included to 
determine the final pumping rate per model cell. This means each DSS well represents three persons. 
Applied to all four counties: 


e Collier County = 123.67 gpd/p x 3 = 371.01 gpd/well 

e Glades County = 114.58 gpd/p x 3 = 343.74 gpd/well 

e Hendry County = 67.43 gpd/p x 3 = 202.29 gpd/well 

е Lee County = 71.05 gpd/p x 3 = 213.15 gpd/well 

From information received from Lee County and discussions with SFWMD Water Supply Planning staff, 
two additional pieces of information were used to refine the DSS pumping rates. First, the DSS well average 
PCUR for each county was increased by a factor of one to account for private landscape irrigation. Applied 
to all four counties: 


Collier County = 371.01 gpd/well + 123.67 = 494.68 gpd/well 
Glades County = 343.74 gpd/well + 114.58 = 458.32 gpd/well 
Hendry County = 202.29 gpd/well + 67.43 = 269.72 gpd/well 
Lee County = 213.15 gpd/well + 71.05 = 284.20 gpd/well 


Second, the pumping rates for all DSS wells were adjusted over time according to seasonal variations 
measured at DSS wells in a portion of Lee County. This information was obtained from a Lehigh Acres 
groundwater flow model being developed by A.D.A. Engineering, Inc. A.D.A. Engineering, Inc. provided 
the SFWMD with calibrated DSS well pumping rates based, in part, on private landscape irrigation use and 
monthly rate variations over the course of a 2-year period (January 2013 through December 2014). The 
monthly variation in DSS well pumping rates was used to derive a multiplier (ranging from 0.86 to 1.14) 
that was applied to the PCUR for each DSS well for each month of the LWCSIM calibration period (January 
1999 through December 2012). 


As a check on the above water extraction rate per well, the overall, four-county average PCUR of 
373 gpd/well divided by 3 persons/well = 124 gpd/p. This number is almost identical to the 2015 finished 
residential water demand value in the WaSUP database for Collier County (123.67 gpd/p). 


The final PCUR per well per county was assigned to the DSS wells in the model cells, with adjustment for 
seasonal variations applied monthly from January 1999 through December 2012. 


Figure C-2 shows the locations of all DSS wells included in the LWCSIM. There are 59,340 DSS wells in 
Lee, Collier, Hendry, and Glades counties, represented by a total of 22,139 model cells. 
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Figure C-2. Domestic Self-Supply well coverage in the model. 


This DSS well coverage and an Excel pivot table were used to prepare the spreadsheet shown in Table C-1, 
which was used to build part of the WEL package for the LWCSIM. Using the procedure described earlier, 
DSS well data, including layer, row, column, average PCUR, and monthly per well PCUR, were entered 
into this spreadsheet. The data were used in the WEL package for transient MODFLOW flow modeling. 
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Table C-1. | Domestic Self-Supply well data. 


won O) AUNE 


re 
но 


А [В [С | р | Е Е | G | H | I | J | K | L | 
Layer ROW COL Q gpdpw Q mgdpw Jan-99 Feb-99 Mar-99 Apr-99 May-99 Jun-99 201-99 
1 110 174 163.4991534 0.000163499 0.000168417 0.000187956 0.000169167 0.000176838 0.000166456 0.00014418 0000141154 
1 110 177 176.5850443 0.000176585 0.000181896 0.000203 0.000182707 0.000190992 0.000179779 0.00015572| 0.000152451 
1 111 178 169.56996 0.00016957| 0.00017467 0.000194935 0.000175449 0.000183405 0.000172637 0.000149534 0.000146395 
1 112| 137 195.965799 0.000195966) 0.00020186 0.000225279 0.00020276 0.000211954 0.00019951 0.000172811 0.000169183 
1 112 171 159.9002048 0.0001599 0.00016471 0.000183819 0.000165444 0.000172946 0.000162792| 0.000141007 0.000138047 
1 112 174 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387, 0.00028934 0.00025062 0.000245359 
|1 113 159 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| (0.00025062 0.000245359 
1 113 160 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934 0.00025062 0.000245359 
1 113 168 146.944195 0.000146944 0.000151364 0.000168925 0.000152039 0.000158933 0.000149602 0.000129582 0.000126861 
1 113 173 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| 0.00025062| 0.000245359 
1 113 187 2842 0.0002842| 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934 0.00025062 0.000245359 
|1 114 159 284.2 0.0002842| 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934 0.00025062 0.000245359 
1 114 164 142.0923451 0.000142092 0.000146366 0.000163347 0.000147019 0.000153685 0.000144662 0.000125303 0.000122673 
1 114 166 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| (0.00025062 0.000245359 
|1 114 169 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934 (0.00025062 0.000245359 
1 114 170 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| 0.00025062 0.000245359 
1 115 167 28.81442074 2.88144Е-05 2.96811Е-05 3.31246Е-05 2.98134E-05 3.11653E-05 2.93355Е-05 2.54098E-05 2.48764Е-05 
1 115 171 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934 0.00025062 0.000245359 
1 115 177 158.4359975 0.000158436 0.000163201 0.000182136 0.000163929 0.000171362 0.000161301 0.000139716 0.000136783 
1 115 186 2842 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| 0.00025062 0.000245359 
1 115 187 159.8681165 0.000159868 0.000164677 0.000183782 0.000165411 0.000172911 0.000162759 0.000140978 0.000138019 
1 116 137 484.8676176 0.000484868 0.000499451 0.000557397 0.000501677 0.000524426 0.000493637 0.000427577 0.000418601 
1 116 159 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| (0.00025062 0.000245359 
1 116 177 146.129343 0.000146129 0.000150525 0.000167988 0.000151195 0.000158051 0.000148772 0.000128863 0.000126158 
1 117 137 600.6779659 0.000600678| 0.000618745] 0.00069053 0.000621503 0.000649685 0.000611541 0.000529703 0.000518584 
1 117 138 206.4800919 0.00020648 0.000212691 0.000237366 0.000213638 0.000223326 0.000210214 0.000182083 0.000178261 
1 117 161 2842 0.0002842 0.000292748 0.000326712 0.000294053, 0.000307387, 0.00028934) 0.00025062 0.000245359 
1 117 162 852.5817742 0.000852582 0.000878226 0.000980115 0.00088214 0.000922141 0.000868001 0.000751843 0.00073606 
1 117 163 284.2 0.0002842 0.000292748 0.000326712 0.000294053| 0.000307387 0.00028934| (0.00025062 0.000245359 
1 117 167 284.2 0.0002842 0.000292748 0.000326712 0.000294053 0.000307387 0.00028934| 0.00025062 0.000245359 
1 117 178 147.0996218 0.0001471 0.000151524 0.000169104 0.000152199 0.000159101 0.00014976 0.000129719 0.000126996 
1 117 180 142.090377 0.00014209| 0.000146364 0.000163345 0.000147016 0.000153683 0.00014466 0.000125301 0.000122671 


Іп summary, the steps used to arrive at the input table (Table C-1) for the MODFLOW WEL package аге 
as follows: 


Gather DSS well locations, casing depths, and total well depths from each county. County health 
departments are a good starting point for this information. Lee and Collier counties were able to 
provide GIS shapefiles and Excel spreadsheets with enough information to locate all DSS wells. 


Create a GIS point coverage of all DSS wells. Intersect the DSS points with the LWCSIM grid to 
assign a row and column to each DSS well. 


Run this coverage and the model layer coverage through a GIS analysis tool to assign layers to each 
DSS well point. 


Develop a PCUR for DSS wells in each county. 


Create a master DSS well Excel spreadsheet with the following fields: row, column, county, X, Y, 
case depth, total depth, PCUR per well, PCUR per layer, top of uppermost layer tapped, bottom of 
lowermost layer tapped. 


Create a pivot table in this spreadsheet that sums the PCURs in cells containing multiple DSS wells 
and assigns a PCUR to each model cell. 


Create a MODFLOW WEL Package spreadsheet that contains the following fields: layer, row, 
column, PCUR (gpd/well), PCUR (mgd/well), and a monthly time series of PCURs for transient 
flow model runs. 


LITERATURE CITED 
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APPENDIX D: 
RIVER STAGE MISSING DATA GAP FILLING 


For this modeling effort, incomplete stage data at multiple monitoring stations were filled using two 
methods: 1) a simple averaging, and 2) a sequence-based Markovian stochastic model (Assegid and Melesse 
2011). 


1. For stations with incomplete stage data of 2 years or less, the monthly average of 1999 to 2014 was 
used to fill the gap. 


2. For stations with incomplete stage data of greater than 2 years, a sequence-based Markovian 
stochastic model (Assegid and Melesse 2011) was used. The method was tested for groundwater 
fluctuation simulation in South Florida and extended to surface water level prediction using the 
same analogy of groundwater fluctuation to the stage level fluctuation simulation. The result was 
presented to hydrologists and a consensus was reached on its use. Accordingly, the model was used 
to fill missing data for three river monitoring stations and three shallow groundwater stations 
(Table D-1). The shallow groundwater stations are in close proximity to ungauged rivers and thus 
serve as proxies for water level in the absence of monitoring. 


Table D-1. Missing data by station. 


Station Status Missing Data Years Duration (years) Root Mean Square Error 
VALI75 Stage 1999 — 2005 6 0.4 

HF3 G Well 1999, 2013 – 2014 1,2 1.2, 1.1 
WFI_G Well 1999, 2003 — 2014 1,12 1.0, 1.1 


The stochastic method simulated monthly water level fluctuations based on available data. The data were 
used to develop and test the stochastic model. The parameters were later used to predict the period of time 
missing data were observed. The model used the least root mean square error and maximum likelihood ratio 
of 90% at the time of simulation and testing. In all cases of the forecast, the variance of the observed was 
preserved to uphold the stationarity assumption of the Markovian approach. The result suggests as an 
acceptable error for all cases: a root mean square error of 0.1 to 1.5 ft (3% to 5% of the head). A similarity 
of patterns between the observed and simulated outputs was observed, as indicated in Figures D-1 to D-3. 
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Figure D-1. Observed versus predicted values for the missing data gap at the WF1_G shallow well. 
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Figure D-2. Observed versus predicted values for the missing data gap at the HF3_G shallow well. 
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Figure D-3. Observed versus predicted values for the missing data gap at the VALI75 gauging station. 
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APPENDIX Е: 
DRAIN COVERAGE PROCEDURE 


E-1 


All shapefiles and feature classes illustrated in this appendix as well as the ArcMap projects used to 
assemble the coverages are stored on the South Florida Water Management District' s (SFWMD’s) server 
at: \\ad.sfwmd.gov\dfsRoot\data\wsd\MOD\RCE. 


Two coverages were needed to prepare the drains coverage for MODFLOW: 1) National Hydrography 
Dataset (http://nhd.usgs.gov/) medium resolution (NHDM), and 2) the digital elevation model (DEM) for 
South Florida, light detection and ranging (LIDAR) 100 feet (ft), found in the SFWMD’s Geospatial Data 
Library. 


The NHDM was clipped to the Lower West Coast Surficial and Intermediate Aquifer Systems Model 
(LWCSIM) domain and saved as a feature class. After eliminating all but StreamRiver, CanalDitch, and 
ArtificialPath from the coverage, the remaining drain flow polylines in the LWCSIM domain are shown in 
Figure E-1. 
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Figure E-1. Drains in the model area. 
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Intersect the drain flow polylines with the LWCSIM 1,000 ft x 1,000 ft grid mesh: 


arcpy.Intersect analysis(in featuresz " NHDMflow clipped 101215 #;'LWCSIM Model Grid 
Grp/LWCSIM mesh' £",out feature class- "/ad.sfwmd.gov/dfsRoot/data/wsd/MOD/RCE/LWCSIM 
DRNs.gdb/NHDMflow clipped intMesh 101215" join attributes- "ALL", cluster. tolerance- "output t 
ype="INPUT") 


Results in flow line segments inside model cells (Figure E-2). 
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Figure E-2. Drains inside model cells. 
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The next coverage used was the DEM for South Florida, LIDAR 100 ft (DEM E(South Florida) Y(2013) 
S(FDEM, SFWMD) (LiDAR, 100 ft)), loaded from the SFWMD’s Geospatial Data Library. This raster 
coverage was clipped to the active model domain boundary (Figure E-3). 
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Figure E-3. Digital elevation model within the model area. 


The clipped DEM was integerized using the Raster Calculator to multiply “NAVD88Raster” by 100. Then 
Extract by Mask was used on the previous two coverages (Figure E-4). 


# The following inputs are layers or table views: "DEM E(South Florida) Y(2013) $(ЕРЕМ, SFWMD) 
(LiDAR, 100 feet)", "NHDMflow clipped intMesh 101215" 


arcpy.gp.ExtractByMask sa("DEM E(South Florida) Y(2013) <(ЕРЕМ, SFWMD) (LiDAR, 100 


feet)", NHDMflow. clipped intMesh 101215", "//ad.sfwmd.gov/dfsRoot/data/wsd/M OD/RCE/LWCSIM 
DRN 02. gdb/FlowLineDEMInt") 
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Figure E-4. Drain segments in model cells with elevation data added. 
The Aggregate tool was used to calculate the mean of the integerized elevations (Figure E-5). 
# The following inputs are layers or table views: "FlowLineDEMInt" 


arcpy.gp.Aggregate_sa("FlowLineDEMInt", "//ad.sfwmd. gov/dfsRoot/data/wsd/MOD/RCE/LWCSIM DRN 
02. gdb/DrnCellDEM","10","MEAN","EXPAND","DATA") 
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Figure E-5. Model cell drain segments with average elevation data. 
The mean elevations were intergerized again using the Raster Calculator (Figure E-6). 


arcpy.gp.RasterCalculator. sa("""Int("DrnCellDEM")""", "//ad.sfwmd. gov/dfsRoot/data/wsd/MOD/RCE/L 
WCSIM DRN 02. gdb/DrnCellDEMint") 
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Figure E-6. Model cell drain segments with average elevation data integerized. 

The Raster to Polygon tool was used to create polygons of the model drain segments (Figure E-7). 

# The following inputs are layers or table views: "DrnCellDEMint" 

arcpy.RasterToPolygon conversion(in rasterz "DrnCellDEMint",out. polygon features- "//ad.sfwmd.go 


v/dfsRoot/data/wsd/MOD/RCE/LWCSIM DRN 
02.gdb/DRNcellDEMPoly",simplify- "NO. SIMPLIFY",raster. field "Value") 
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Figure E-7. Model drain cell segments converted to polygons. 


The Calculate Field tool was used on this coverage to create a new field in the attribute table called 
DRN_Elev. 


# The following inputs are layers or table views: "DRNcellDEMPoly" 


arcpy.CalculateField_management(in_table="DRNcellDEMPoly" field="DRN_Elev", expression="[grid 
code ]/100", expression_type="VB",code_block="#") 


The previous steps were important because the elevation data from the DEM raster were needed to assign 
elevations to each drain flow segment. 
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The coverage NHDMflow_clipped_intMesh_101215 was intersected with DRNcellDEMPoly to create 
FlowLineElevRoCo, which contains the drain flow lines inside mesh cells, drain elevations, and 
row/column (Figure E-8). 


arcpy.Intersect_analysis(in_features="NHDMflow_clipped_intMesh_101215 #; DRNcellDEMPoly 
#" out_feature_class="//ad.sfwmd. gov/dfsRoot/data/wsd/MOD/RCE/LWCSIM DRN 
02. gdb/FlowLineElevRoCo" join attributes- "ALL",cluster. tolerancez "f", output type- "INPUT") 
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Figure E-8. Final drain segments with elevation data and row/column. 


E-9 


This was intersected with LWCSIM_Soils_082515 to create the final coverage (Figure E-9), which contains 
all the pertinent data for the DRN package (row, column, segment length, drain elevation, and soil texture). 
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Figure E-9. Final drain segments with soil textures added. 


arcpy.Intersect_analysis(in_features="DRN_FlowSumLenRoCo_MeanElev #; ‘Soils 
Grp/LWCSIM_Soils_082515' 

#" out_feature_class="//ad.sfwmd. gov/dfsRoot/data/wsd/MOD/RCE/LWCSIM DRN 

02. gdb/Final_DRN_Coverage_101415" join_attributes="ALL",cluster_tolerance="#",output_type="INP 
UT") 
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The River coverage was subtracted from the final drain coverage to eliminate overlap with stream features 
that are considered Rivers in the LWCSIM (Figure E-10). 
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Figure E-10. Drains without rivers. 


The features were exported to create the Excel spreadsheet shown in Table E-1, which was used to generate 
the DRN package for MODFLOW. 


Table E-1. Drains coverage information. 


A B C D E F G H I J K L 
LAYER ROW COL ROWCO ELEV f msl WIDTH f THICKf © SEGLENGTHf CONDUCTANCEf//d Texture K f/d 

! H 548 425 7 548425 0.1 30 1 255.7869372 107430.5136 FS 14 
! 1 548 426 " 548426 1.48 30 1 1017.369013 427294.9855 FS 14 

1 548 427 " 548427 1.48 30 1 501.5893038 210667.5076 FS 14 
' 1 548 430 " 548430 1.76 30 1 1029.448798 432368.4952 FS 14 
i 1 547 421 ” 547421 -0.76 30 1 416.9044981 175099.8892 FS 14 
К 1 547 422 " 547422 -0.76 30 1 36.59626064. 15370.42947 FS 14 
i 1 547 423 ” 547423 0.35 30 1 413.9235281 173847.8818 FS 14 
! 1 547 424 " 547424 0.08 30 A 339.6062004 142634.6042 FS 14 
D 1 547 427 ” 547427 1.48 30 1 647.7809808 272068.0119 FS 14 
1 1 547 428 ” 547428 1.18 30 1 1108.805493 465698.3071 FS 14 
2 1 547 429 7 547429 1.18 30 1 306.609305 128775.9081 FS 14 
3 1 547 430 " 547430 1.18 30 1 999.4502819 419769.1184 FS 14 
4 1 546 421 " 546421 0.81 30 1 983.9156952 413244.592 FS 14 
5 1 546 422 " 546422 0.81 30 1 1434.51678 602497.0476 FS 14 
5 1 546 423 ” 546423 0.51 30 1 1098.716915 461461.1043 FS 14 
7 1 546 425 " 546425 0.47 30 1 367.1718938 154212.1954 FS 14 
8 1 546 429 ” 546429 1.12 30 1 1295.863012 544262.465 FS 14 
9 1 546 430 " 546430 1.12 30 1 147.2199574 61832.38211 FS 14 
2 1 545 421 " 545421 0.03 30 1 978.7979272 411095.1294 FS 14 
1 1 545 422 7 545422 0.3 30 1 100.4156987 42174.59345 FS 14 
2 1 545 423 " 545423 0.33 30 1 323.1706273 135731.6635 FS 14 
3 1 545 424 " 545424 0.33 30 1 1993.125476 837112.6999 FS 14 
4 1 545 425 7 545425 0.33 30 1 2159.429875 906960.5475 FS 14 
5 1 545 426 ” 545426 0.79 30 1 1240.823006 521145.6625 FS 14 
5 1 545 427 " 545427 1.1 30 1 1097.421579 460917.0632 FS 14 
7 1 545 430 " 545430 1.22 30 1 1707.336569 717081.359 FS 14 


The LAYER, WIDTH f, and THICK f columns were added in Excel. All flow lines are in layer 1. Width 
and thickness for all flow line segments were estimated to be 30 ft wide and 1 ft thick of drain bed material. 
These data, along with soil textures, were used to estimate hydraulic conductivity and conductance. 
Hydraulic conductivity was based on soil texture and derived from literature. Conductance was calculated 
as: 


C=KLW/M 


Where: 


C is the conductance (ft?/day) 

K is hydraulic conductivity (ft/day) 

L is length of the drain segment (ft) 

W is width of the drain segment (estimated to be an average of 30 ft) 

M is thickness of the drain bed sediments (estimated to be an average of 1 ft) 


A final check on drain elevations was completed by comparing Marco Water Engineering, Inc. (2006) weir 
crest elevations with corresponding topographic elevations (Figure E-11). 
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Figure E-11. Model drains and Marco Water Engineering, Inc. (2006) weir elevation data. 


The average difference in elevation between land surface and the weir crest elevations was 3.4 ft lower. 
Therefore, an adjustment of 4 ft was applied to all drain segments in the final drain coverage. 
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APPENDIX F: 
WETLAND COVERAGE PROCEDURE 
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This appendix describes the geographic information systems (GIS) procedure used to produce the WTL 
package from 2010 land use/land cover (LULC) data. LULC years 2000, 2004, and 2014 were used in 
similar procedures to produce WTL packages for other model years. All shapefiles and feature classes 
illustrated in this appendix as well as the ArcMap projects used to assemble the coverages are stored on the 
South Florida Water Management District’s (SFWMD’s) server at: 
\\ad.sfwmd.gov\dfsRoot\data\wsd\MOD\RCE. 


For Wetland Year 2010: 


1. In the 2010 WTL Coverage Grp, begin with LWCSIM_Wetlands_2010. This is a feature class 
containing LULC) data that have been reduced to only wetlands categories. These have been broken 
out to Levels 1, 2, and 3 and clipped to the model domain (purple boundary in Figure F-1). The 
red box is a closeup area used for the remaining figures in this appendix. 
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Figure F-1. Overview of wetland coverage. 


2. Closeup area (Figure F-2): 
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Figure Е-2. Closeup. 
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3. Dissolve to the field LEVOI, no multipart features. This step creates contiguous wetland polygons 
within the model area. Save the output as wetl2010 І.ЕУ0І Diss (Figure F-3). 


- - Em 
— „А е Iw Ж CIEL ee TA 
ПӘ $E! =| Ё e ao 1 в 
= Layers у 
& L1 2000 WTL Coverage Grp 
& O 2004 WTL Coverage Grp 
В M 2010 WTL Coverage Grp 
ш О GrdSelO2IntWTLO2 
в О GrdSel02 
ш O GrdSel01Diss 
в O GrdSel01 
я O wetl2010 disLev02 IntMesh 500K 
ш O мей2010 disLev02 IntMesh 
& О wetl2010 disLev02 
в O wetl2010grdcls dis 3m intLU 
& O wetl2010grdcls dis 3m 
в O wetl2010grdcls dis 
& O wetl2010grdcls 
& O wetl2010LVO1DisMesh Diss 500K 
& О wetl2010LVO1DisMesh Diss 
mu Weti2010LVO1Dis mesh 


| 
же 
Е 
| 


& O LWCSIM Wetlands 2010 
& [1 Marco Wetland Coverage 
ш O FLNWI Grp 
в М 055 Wells Grp 
& М LWCSIM Model Grid Grp 
& М SFWMD Basemap Grp 
& O LWCSIM AHED AII Canals clip 
& LJ MedStdDev WTL Pts 
& [1 Apps2010DivLivAgrNur 


Figure F-3. Contiguous wetland polygons. 
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4. Take wetl2010_LEV01_Diss and intersect it with the Lower West Coast Surficial and Intermediate 
Aquifer Systems Model (LWCSIM) grid (LWCSIM_mesh in LWCSIM Model Grid Grp). The 
result is a set of polygons with model grid row and column attributes for any wetland polygons in 
the model area. Save the output file as wetl2010LV01Dis_mesh (Figure F-4). 
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Figure F-4. Contiguous wetland polygons intersected with model grid. 
5. Take wetl2010LVO1Dis mesh and dissolve to field FID LWCSIM mesh (this is the feature ID # 


assigned to the model grid cells). “Yes’ to multipart, dissolve lines. Model cells with multiple 
wetland polygons are merged. This file should be saved as wetl2010LVO01DisMesh Diss. 
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6. Take wetl2010LV01DisMesh_ Diss and select by attribute Shape area > 500,000. This identifies 
model cells that are more than 50% wetland (500,000 is 50% of 1,000 square foot cells). Export 
the data to create wetl2010LVO1DisMesh, Diss 500K (Figure Е-5). 
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Figure F-5. Model cells containing 50% ог more wetlands. 
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7. From LWCSIM_mesh, select by location against the coverage created in step 5 
(wetl2010LV01DisMesh_Diss_500K) as the source. Spatial selection method for the target layer 
feature(s): contain (Clementini) the source layer features. This will select only the cells in 
LWCSIM mesh that contain (or touch) the wetlands polygons created in step 5. Export the data to 
the file wetl2010grdcls. This coverage contains only the cells with wetland data for the model 
(Figure F-6). 
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Figure F-6. Grid cells representing wetlands. 


8. Dissolve wetl2010grdcls to create wetl2010grdcls dis. This will join together adjacent grid cells to 
form continuous polygons (Figure F-7). 
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Figure F-7. Wetland grid cells dissolved into contiguous polygons. 


9. From wetl2010grdcls dis, select by attribute the polygons that are > 3 million (3 grid cells in area), 
and export the data to wetl2010grdcls_dis_3m (Figure F-8). 
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Figure F-8. Wetland grid cell polygons larger than three grid cells in size. 


10. Intersect wetl2010grdcls_dis_3m with the original LWCSIM_Wetlands_2010 coverage to create 
wetl2010grdcls_dis_3m_intLU. This restores the original wetlands data to only the data within the 
model grid polygon areas arrived at in step 8 (Figure F-9). 
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Figure F-9. Original wetland data restored to wetland grid polygons. 
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11. Take wetl2010grdcls_dis_3m_intLU and dissolve to the field LEVO3. Output this to the file 
wetl2010_disLev03. This will reduce the wetlands data to Level 3 wetland categories. 


F-10 


12. Take wetl2010_disLev03 and intersect with the model grid (LWCSIM_mesh). Send the output to 
the file wetl2010_disLev03_IntMesh. This now shows only the portions of Level 3 wetlands that 
intersect model grid cells (Figure F-10). 
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Figure F-10. Model grid cells containing only Level 3 wetlands. 
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13. Take wetl2010 disLev0O3 IntMesh, select by attribute shape area > 500,000 and export the data to 
create the file wetl2010 disLev03 IntMesh 500K. This step will eliminate all wetlands polygons 
smaller than 50% of a model grid cell. However, individual groups of less than three cells remain. 
The final steps will eliminate those and leave only groups of three or more cells, each of which 
contains 50% or more wetland area (Figure F-11). 
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Figure F-11. Model grid cells containing 5096 or more Level 3 wetlands. 
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14. From the mesh, select by attribute the grid cells that match with information in step 12 to create 
GrdSel03 (Figure F-12). 
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Figure F-12. Resulting selection with model grid. 
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15. Dissolve the cells to create contiguous polygons and save to GrdSel03Diss (Figure Е-13). 


[1 2000 WTL Coverage Grp 
[1 2004 WTL Coverage Grp 
М 2010 WTL Coverage Grp 
I О GrdSelO4IntWTLO4 

E О GrdSelO3IntWTLO3 

E О GrdSelO3Diss 3m 
EM 


m 
Æ О GrdSel03 
E ІП GrdSelO2IntWTLO2 
Ею О GrdSel02 
ЕП GrdSel01Diss 
Е О GrdSel01 
Ею О wetl2010_disLev03_IntMesh_500K 
ЕП wetl2010 disLev02 IntMesh 500K 
E О wetl2010 disLev03 IntMesh 
E О wetl2010 disLev02 IntMesh 
E O wetl2010 disLev03 
# О wetl2010 disLev02 
E О wetl2010grdcls dis 3m intLU 
Е О wetl2010grdcls dis 3m 
Е О wetl2010grdcls dis 
Е О wetl2010grdcls 
Е О wetl2010LVO1DisMesh Diss 500K 
I О wetl2010LVO1DisMesh Diss 
= О wetl2010LVO1Dis mesh 
Б О wetl2010 LEVO1 Diss 
Е О LWCSIM Wetlands 2010 
O Marco Wetland Coverage 


m 
D FLNWI Grp 
М DSS Wells Grp 
li ІП DSS Master 081715 for SHP 
Is ІП DSS Master 081415 for SHP 


ЕП ColloierCounty address 


Æ О LeeCountyPermittedWells2014022 


г | 


Collier 


Figure F-13. Model grid selection dissolved to contiguous polygons. 
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16. Select only the polygons that are three cells in size (3 million square feet) or more to create 
GrdSel03Diss_3m (Figure F-14). 
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Figure F-14. Resulting grid polygons reduced to > grid cells in size. 
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17. Finally, intersect the polygons with the original wetlands LULC coverage to create 
GrdSel04IntWTL04. The attribute table of GrdSel04IntWTL04 will contain the critical wetland 
data for each active cell in the model grid (Level 3 wetland category, row, and column) 
(Figure F-15). 
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Figure F-15. Intersection with original wetland data. 


This coverage was intersected with soils and model layers to produce the Excel spreadsheet shown in 
Table F-1. 


Table F-1. | Wetland coverage information. 


А вос D Е F G H 1 Й к L M 
ROW COL Layer ROWCO Topo NGVD2  FLUCCSLEVO3 FLUCCS Description “REG CAT  MaxSoilArea Percent Cell Soil MUNAME TEXTURE  KValuef/d _ 

362 490 1762490 9.74336 6410 Freshwater marshes 2 3154536724 0.32 Riviera-Copeland-Boca (s1593) SAND 141 
442 29 1742279 1.0626 6120 Mangrove swamps 1 17936.72304 1.79 Water-Terra Ceia-Perrine-Pennsuco-Okeelanta (51582) MUCK 0.28 
363 490 1763490 9.73396 6410 Freshwater marshes 2 43613.99505 4.36 Riviera-Copeland-Boca (51593) SAND 141 
365 490 1765490 9.71468 6410 Freshwater marshes 2 98752.57081 9.88 Riviera-Copeland-Boca (51593) SAND 141 
364 490 1764490 9.72752 6410 Freshwater marshes 2 119836.8069 11.98 Riviera-Copeland-Boca (51593) SAND 141 
285 450 1785450 15.4533 6210 Cypress 2 133869.7747 13.39 RIVIERA SAND, LIMESTONE SUBSTRATUM, DEPRESSIONAL SAND 141 
257 394 1757394 22.9003 6170 Mixed wetland hardwoods 2 135266.9386 13.53 HOLOPAW SAND, DEPRESSIONAL SAND 141 
257 391 1757391 23.1885 6210 Cypress 2 147648.7544 14,76 CHOBEE FINE SANDY LOAM, LIMESTONE SUBSTRATUM, DEPRESSIONAL FSL 8 
275 391 1575391 20.0062 6240 Cypress-pine-cabbage palm 3 154449.4619 15.44 CHOBEE FINE SANDY LOAM, LIMESTONE SUBSTRATUM, DEPRESSIONAL FSL 8 
267 405 1767405 19.9912 6170 Mixed wetland hardwoods 2 1611767175 16.12 JUPITER-OCHOPEE-ROCK OUTCROP COMPLEX в 14 

17 204 1017204 43.1661 6300 Wetland forested mixed - 162978.6515 16.30 MALABAR FINE SAND, HIGH FS 14 
233 189 1733189 16.3137 6210 Cypress 2 164669.2076 16.47 PINEDA FINE SAND FS 14 
230 424 1730424 24.0095 6410 Freshwater marshes 2 1660865379 16.61 BASINGER SAND. SAND 141 
372 249 17372249 742611 6200 Wetland coniferous forests - 166380.2917 16.64 OLDSMAR FINE SAND, LIMESTONE SUBSTRATUM FS 14 
183 310 1783310 28.8782 6170 Mixed wetland hardwoods 2 167204.3198 16.72 СНОВЕЕ FINE SANDY LOAM, DEPRESSIONAL ғы в 
307 272 17307272 14.2904 6250 Hydric pine flatwoods 3 1701739114 17.02 HOLOPAW FINE SAND FS 14 
260 199 1760199 15.7556 6250 Hydric pine flatwoods 3 170289.9024 17.03 PINEDA FINE SAND, DEPRESSIONAL в 14 
322 478 1722478 12.7938 6170 Mixed wetland hardwoods 2 1721444775 17.21 JUPITER FINE SAND FS 14 
322 478 1522478 12.7938 6170 Mixed wetland hardwoods 2) 1721444775 17.21 Riviera-Copeland-Boca (51593) SAND 141 
219 244 1719244 29.4802 6210 Cypress 2 1731348112 17.31 FELDA FINE SAND, DEPRESSIONAL FS 14 
295 485 1795485 12.7506 6410 Freshwater marshes 2 174598.5322 17.46 CHOBEE SANDY LOAM, DEPRESSIONAL я 2 
212 452 1512452 16.8457 6210 Cypress 2 175447.9349 17.54 HOLOPAW SAND SAND 141 
448 272 1740272 3.18447 6120 Mangrove swamps 1 177454.4447 17.75 DURBIN AND WULFERT MUCKS, FREQUENTLY FLOODED миск 0.28 
186 405 1786405 27.0024 6430 Wet prairies 3 179236.1882 17.92 WINDER FINE SAND FS 14 
218 199 1718199 20.7428 6300 Wetland forested mixed - 179426.6203 17.94 VALKARIA FINE SAND, DEPRESSIONAL FS 14 
293 329 17793329 16.7803 6170 Mixed wetland hardwoods 2 180030.4873 18.00 CHOBEE, WINDER, AND GATOR SOILS, DEPRESSIONAL FSL 8 
296 309 1796309 15.594 6210 Cypress 2 1802949219 18.03 BOCA FINE SAND в 14 
447 22 1047272 2.98785 6120 Mangrove swamps 1 1809234907 18.09 DURBIN AND WULFERT MUCKS, FREQUENTLY FLOODED MUCK 0.28 
291 337 17791337 16.4267 6300 Wetland forested mixed - 183232.8896 18.32 HILOLO, JUPITER, AND MARGATE FINE SANDS FS 14 
205 264 1205264 25.4073 6410 Freshwater marshes 2. 183596648 18.36 TUSCAWILLA FINE SAND в 14 
331 325 1331325 15.9546 6250 Hydric pine flatwoods 3 1842347491 18.42 OCHOPEE FINE SANDY LOAM FSL 8 
оқа ABS 4 Зколас 17 аяз BAIN Erachwatar marchac I 190262 £124 ая ла DI ANTATION миск миск non 


Data from the spreadsheet were used to generate the array files needed for the WTL package for 
MODFLOW. 
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APPENDIX С: 
LANDSCAPE IRRIGATION DATA DEVELOPMENT PROCEDURE 


Identify landscape irrigation (LSI) locations and calculate flow rates 


a. Delineate developed areas in geographical information systems (GIS) land use coverages. 


b. Overlay utility service area (potable and reclaimed water) map and identify different 
application regions belonging to each county. 


c. Overlay wetland/water features coverage and percent impervious coverage to delineate 
pervious areas within each cell and each irrigation application area. 


Calculate monthly application rates. 


a. Using the map produced in step 1, calculate the percent pervious areas in each cell located in 
developed areas in each irrigation application area in a given county. 


b. Normalize the pervious area calculated for each cell by the total pervious area within the 
respective irrigation application area in a given county. 


c. Multiply the total irrigation given by the normalized percent pervious area for each cell to 
calculate the flow rate for each cell. This is the average annual daily flow; therefore, the 
calculated LSI is the average annual daily LSI flow. 


хв Ai xfi 


num_class 


LSI rate in cell = Group LSI rate * ———_______ 
All cells in group Di Aj * fj 


Develop the temporal variation curve for LSI. 


a. Peaking factors were derived from the outdoor portion of the average pumping of all utilities 
in a given county (flow-weighted average). Utilities with larger pumping volume have more 
influence on the peaking factors than utilities with smaller pumping volume. The following 
steps describe the derivation of peaking factors and temporal variation curve for LSI. 


i. Locate utilities that provide service within a given county. 


ii. Calculate the total pumping volume for each of the stress periods from all wells/utilities 
for each county. 


iii. Obtain the outdoor portion of the total pumping volume used for LSI for each stress period 
by subtracting the indoor water use from values obtained in step ii. 


iv. The array of peaking factors for each county is obtained by dividing the outdoor potable 
LSI county array obtained in step iii by the total outdoor portion of the water use. 


v. Calculate the temporal variation of LSI in each cell by multiplying the average annual daily 
flow by the peaking factors obtain in step iv. 


vi. Sum all the LSI in each cell for the total number of stress periods and calculate the average 
to check the correct distribution of LSI within the irrigable area in each county. 
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4. Identify and calculate the daily irrigation rates. 


a. 


During agricultural demand estimation, LSI demands were calculated using the Agricultural 
Field-Scale Irrigation Requirements Simulation (AFSIRS) demand model. 


Using AFSIRS-generated irrigation demand curves, identify the days irrigation is required. 


For a new cell that was not included in the previous approach, an average irrigation demand 
curve is used for the county where the cell is located. 


Daily application rates are determined by multiplying the monthly average by the number of 
days in a given month and dividing by the number of days that actual irrigation occurs 
(calculated with AFSIRS). Then irrigation is applied on the AFSIRS-calculated days in a given 
month. 


In the model, irrigation is applied to each cell as an irrigation depth, which is obtained by 
dividing the irrigation rates by the model cell area (1,000 x 1,000 square feet). 


Daily irrigation depths are added to the total rainfall array that is run through the Curve Number 
method to partition rainfall plus irrigation into runoff and infiltration. 


The infiltrating portion is sent through the ET-Recharge-Runoff program to develop 
MODFLOW evapotranspiration and recharge data sets. 


The LSI development flow chart is shown in Figure G-1. 
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Input to LSI 


County LSI totals derived from PWS and WWF data 
and GIS coverage 


(51 Algorithm 


Use Input data and calculate monthly 
application rates (MGD) for each model cell 


M2D Pre-Processor 


Develop daily irrigation (in/day) values using monthly application rates from the 151 
algorithm and the irrigation requirement from AFSRIS based calculations 


RAIN Pre-Processor 


Aggregate the daily rainfall values from natural precipitation, AG return flow, and 151 
return flow and create total rainfall (in/day) 


ET/Recharge Algorithm 


Read total Rainfall data (rain*AG return+ LSI ) and 
partitions into runoff (lake, stream, etc.) and 
infiltration (in/day) to unsaturated zone 


MODFLOW ET and Rechargepackage 


Figure G-1. Landscape irrigation data development flow chart. 
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APPENDIX Н: 
FINAL ADJUSTED MODEL PROPERTIES 


The end result of manual transient calibration provided the following property distributions in the model 
layers. In Figures H-1 to H-5, field measurements of horizontal hydraulic conductivity were overlaid to 
show that the calibrated distributions were within an order of magnitude of field measured values 
throughout the model domain. Figures H-6 to H-10 present the transmissivity distribution in each of the 
five aquifer layers. Vertical leakance was adjusted manually for layers 2, 4, 6, and 8, representing the 
vertical leakance between the five aquifer layers. The result of manual adjustment of vertical leakance is 
presented Figures H-11 to H-14. The storage coefficient, or storativity, in the five aquifer layers resulting 
from manual calibration is presented in Figures H-15 to H-19. 
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Figure H-1. Water Table aquifer (layer 1) calibrated hydraulic conductivity. 
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Figure H-2. Lower Tamiami aquifer (layer 3) calibrated hydraulic conductivity. 
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Figure H-3. Sandstone aquifer clastic zone (layer 5) calibrated hydraulic conductivity. 
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Figure H-4. Sandstone aquifer carbonate zone (layer 7) calibrated hydraulic conductivity. 
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Figure Н-5.  Mid-Hawthorn aquifer (layer 9) calibrated hydraulic conductivity. 
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Figure H-6. Water Table aquifer (layer 1) transmissivity from the calibrated model. 
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Figure H-7. Lower Tamiami aquifer (layer 3) transmissivity from the calibrated model. 
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Figure H-8. Sandstone aquifer — clastic zone (layer 5) transmissivity from the calibrated model. 
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Figure H-9. Sandstone aquifer — carbonate zone (layer 7) transmissivity from the calibrated model. 
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Figure H-10. Mid-Hawthorn aquifer (layer 9) transmissivity from the calibrated model. 
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Figure H-11. Calibrated leakance in the Tamiami confining unit (layer 2). 
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Figure H-12. Calibrated leakance in the Upper Hawthorn confining unit (layer 4). 
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Figure H-13. Calibrated leakance in the Sandstone confining unit (layer 6). 
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Figure H-14. Calibrated leakance in the Mid-Hawthorn confining unit (layer 8). 
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Figure H-15. Water Table aquifer (layer 1) calibrated storativity. 
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Figure H-16. Lower Tamiami aquifer (layer 3) calibrated storativity. 
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Figure H-17. Sandstone aquifer-clastic zone (layer 5) calibrated storativity. 
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Figure H-18. Sandstone aquifer-carbonate zone (layer 7) storativity. 
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Figure H-19. Mid-Hawthorn aquifer (layer 9) calibrated storativity. 
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APPENDIX J: 
R-SQUARED VERSUS TNSE FOR LWCSIM CALIBRATION TARGETS 
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The transformed Nash-Sutcliffe efficiency (TNSE) coefficient, a derivative of the Nash-Sutcliffe efficiency 
(NSE) coefficient, was derived for all aquifer layer transient groundwater calibration targets as an additional 
indicator of transient model fit. Through statistical comparison with the R? value, it was determined that a 
TNSE value of 0 roughly corresponds to an К? value of 0.4, allowing for determination of a TNSE 
calibration criterion for the Lower West Coast Surficial and Intermediate Aquifer Systems Model 
(LWCSIM) that is similar to R2. The NSE and TNSE were determined using the following equations: 


XA: — 502 

М$Е=1—-———— 
п 1(0,- 0)? 
і-і1 1 


Where: 


О; = the 1% observed head value at a particular observation well 
S; = the i^ simulated head value at the same well location 


O = the mean of the observed head values at the observation well 


аби 

TNSE = 1 — = 
С) 
і-1 l 


Where: 


Xi = the i" transformed observed head value at a particular observation well 
Xi-0;-0 


X = the mean of the transformed simulated head values at the observation well (this should be zero) 


7; = the i? transformed simulated head value at the same well location 
Zi25-S 

5 = the mean of the simulated head values at the observation well 

5; = the average simulated head at the monitor well 


Oi, 0, and 5; as previously defined 


Nominally, a TNSE of zero corresponds with an R? of 0.4 in that the quantiles for these two statistics are 
about equal at these two thresholds. To be more precise, the R? quantile for an R? of 0.4 is 0.262 (i.e., about 
26% of R? values are less than 0.4). The TNSE quantile for a TNSE of 0.063 is 0.262 (i.e., about 2696 of 
TNSE values are less than 0.063). The TNSE quantile for a TNSE of 0 is 0.252. Based on the above, there 
are three categories of model fit based on TNSE: 


е TNSE < 0: Model quality is less than null-model of observed mean 
e 0<ТМ65Е < 0.5: Low-quality fit 
• 0.5 < TNSE: Moderate- to high-quality fit 


Based on the above, the TNSE was assigned a cutoff of zero. 


The summary of calibration statistics, with the TNSE values included, is presented in Table J-1. 
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Table 1-1. Summary of calibration statistics. 


Calibration 
ee WTA LTA SSA-clastic | ЭЗА- MHA 
carbonate 

Residual mean (<1) -0.34 0.51 0.70 0.98 -0.85 
Error standard deviation (<5) 2.27 2.97 3.22 1.80 3:21 
596 of observed range* 3.29 2.51 3.68 3.26 6.47 
10% of observed range* 6.58 5.02 7.36 6.52 12.94 
Mean absolute error (MAE) (<5) 1.76 2.44 2.38 1.68 2:55 
Root mean square error (<5) 2.29 2,09 3.24 2.01 3,25 
Minimum residual -9,12 -4.12 -9.91 -3.28 -7.38 
Maximum residual 6.45 9,35 5.27 3.78 4,39 
Number of observation points 297 72 29 18 25 
Percentage with MAE < 2.5 ft (50%) 73 61 62 72 56 
Percentage with MAE < 5.0 ft (80%) 97 96 90 100 88 
Percentage with R? > 0.4 (60%) 75 81 76 83 64 
Percentage with TNSE > 0 (60%) 72 85 83 94 64 


LTA = Lower Tamiami aquifer; MHA = Mid-Hawthorn aquifer; R? = Pearson coefficient of correlation; SSA = Sandstone 
aquifer; TNSE = transformed Nash-Sutcliffe efficiency; WTA = Water Table aquifer. 

Green font indicates compliance with calibration criteria. 

Calibration period: 1999 to 2012. 

*Observed head ranges are for information only. 


Figures J-1 to J-5 present the areal distribution of TNSE and R? values for each aquifer layer in the 
LWCSIM. 
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Figure J-1. Areal distribution of transformed Nash-Sutcliffe efficiency coefficient values (left) and В? values (right) for model layer 1. 
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Figure 7-2. Areal distribution of transformed Nash-Sutcliffe efficiency coefficient values (left) and В? values (right) for model layer 3. 
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Figure J-3. Areal distribution of transformed Nash-Sutcliffe efficiency coefficient values (left) and В? values (right) for model layer 5. 
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Figure J-4. Areal distribution of transformed Nash-Sutcliffe efficiency coefficient values (left) and В? values (right) for model layer 7. 


Figure J-5. Areal distribution of transformed Nash-Sutcliffe efficiency coefficient values (left) and В? values (right) for model layer 9. 
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